From f2d061c9210b1495e9e420d7e4cab1dacc9c7adb Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Wed, 24 Jul 2019 22:42:36 -0700 Subject: [PATCH] Test poynting planes on both half-steps --- fdfd_tools/test_fdtd.py | 30 ++++++++++++++++++++++++------ 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/fdfd_tools/test_fdtd.py b/fdfd_tools/test_fdtd.py index 879d3d4..15c5b4d 100644 --- a/fdfd_tools/test_fdtd.py +++ b/fdfd_tools/test_fdtd.py @@ -88,18 +88,19 @@ class BasicTests(): dxes = self.dxes if self.dxes is not None else tuple(tuple(numpy.ones(s) for s in self.epsilon.shape[1:]) for _ in range(2)) dV = numpy.prod(numpy.meshgrid(*dxes[0], indexing='ij'), axis=0) + mx = numpy.roll(self.src_mask, (-1, -1), axis=(0, 1)) + my = numpy.roll(self.src_mask, -1, axis=2) + mz = numpy.roll(self.src_mask, (+1, -1), axis=(0, 3)) + px = numpy.roll(self.src_mask, -1, axis=0) + py = self.src_mask.copy() + pz = numpy.roll(self.src_mask, +1, axis=0) + u_eprev = None for ii in range(1, 8): with self.subTest(i=ii): u_hstep = fdtd.energy_hstep(e0=self.es[ii-1], h1=self.hs[ii], e2=self.es[ii], **args) u_estep = fdtd.energy_estep(h0=self.hs[ii], e1=self.es[ii], h2=self.hs[ii + 1], **args) - mx = numpy.roll(self.src_mask, (-1, -1), axis=(0, 1)) - my = numpy.roll(self.src_mask, -1, axis=2) - mz = numpy.roll(self.src_mask, (+1, -1), axis=(0, 3)) - px = numpy.roll(self.src_mask, -1, axis=0) - py = self.src_mask.copy() - pz = numpy.roll(self.src_mask, +1, axis=0) s_h2e = -fdtd.poynting(e=self.es[ii], h=self.hs[ii]) * self.dt s_h2e[0] *= dxes[0][1][None, :, None] * dxes[0][2][None, None, :] s_h2e[1] *= dxes[0][0][:, None, None] * dxes[0][2][None, None, :] @@ -110,6 +111,23 @@ class BasicTests(): self.assertTrue(numpy.allclose(sum(planes), (u_estep - u_hstep)[self.src_mask[1]]), msg='planes: {} (sum: {})\n du:\n {}'.format(planes, sum(planes), (u_estep - u_hstep)[self.src_mask[1]])) + if u_eprev is None: + u_eprev = u_estep + continue + + s_e2h = -fdtd.poynting(e=self.es[ii - 1], h=self.hs[ii]) * self.dt + s_e2h[0] *= dxes[0][1][None, :, None] * dxes[0][2][None, None, :] + s_e2h[1] *= dxes[0][0][:, None, None] * dxes[0][2][None, None, :] + s_e2h[2] *= dxes[0][0][:, None, None] * dxes[0][1][None, :, None] + planes = [s_e2h[px].sum(), -s_e2h[mx].sum(), + s_e2h[py].sum(), -s_e2h[my].sum(), + s_e2h[pz].sum(), -s_e2h[mz].sum()] + self.assertTrue(numpy.allclose(sum(planes), (u_hstep - u_eprev)[self.src_mask[1]]), + msg='planes: {} (sum: {})\n du:\n {}'.format(planes, sum(planes), (u_hstep - u_eprev)[self.src_mask[1]])) + + # previous half-step + u_eprev = u_estep + class Basic2DNoDXOnlyVacuum(unittest.TestCase, BasicTests): def setUp(self):