diff --git a/fdfd_tools/test/test_fdtd.py b/fdfd_tools/test/test_fdtd.py new file mode 100644 index 0000000..53fe9f8 --- /dev/null +++ b/fdfd_tools/test/test_fdtd.py @@ -0,0 +1,308 @@ +import numpy +import pytest +import dataclasses +from typing import List, Tuple +from numpy.testing import assert_allclose, assert_array_equal + +from fdfd_tools import fdtd + + +prng = numpy.random.RandomState(12345) + +def assert_fields_close(a, b, *args, **kwargs): + numpy.testing.assert_allclose(a, b, verbose=False, err_msg='Fields did not match:\n{}\n{}'.format(numpy.rollaxis(a, -1), + numpy.rollaxis(b, -1)), *args, **kwargs) + +def assert_close(a, b, *args, **kwargs): + numpy.testing.assert_allclose(a, b, *args, **kwargs) + + +def test_initial_fields(sim): + # Make sure initial fields didn't change + e0 = sim.es[0] + h0 = sim.hs[0] + j0 = sim.js[0] + mask = (j0 != 0) + + assert_fields_close(e0[mask], j0[mask] / sim.epsilon[mask]) + assert not e0[~mask].any() + assert not h0.any() + + +def test_initial_energy(sim): + """ + Assumes fields start at 0 before J0 is added + """ + j0 = sim.js[0] + e0 = sim.es[0] + h0 = sim.hs[0] + h1 = sim.hs[1] + mask = (j0 != 0) + dV = numpy.prod(numpy.meshgrid(*sim.dxes[0], indexing='ij'), axis=0) + u0 = (j0 * j0.conj() / sim.epsilon * dV).sum(axis=0) + args = {'dxes': sim.dxes, + 'epsilon': sim.epsilon} + + # Make sure initial energy and E dot J are correct + energy0 = fdtd.energy_estep(h0=h0, e1=e0, h2=h1, **args) + e0_dot_j0 = fdtd.delta_energy_j(j0=j0, e1=e0, dxes=sim.dxes) + assert_fields_close(energy0, u0) + assert_fields_close(e0_dot_j0, u0) + + +def test_energy_conservation(sim): + """ + Assumes fields start at 0 before J0 is added + """ + e0 = sim.es[0] + j0 = sim.js[0] + u = fdtd.delta_energy_j(j0=j0, e1=e0, dxes=sim.dxes).sum() + args = {'dxes': sim.dxes, + 'epsilon': sim.epsilon} + + for ii in range(1, 8): + u_hstep = fdtd.energy_hstep(e0=sim.es[ii-1], h1=sim.hs[ii], e2=sim.es[ii], **args) + u_estep = fdtd.energy_estep(h0=sim.hs[ii], e1=sim.es[ii], h2=sim.hs[ii + 1], **args) + delta_j_A = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii-1], dxes=sim.dxes) + delta_j_B = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii], dxes=sim.dxes) + + u += delta_j_A.sum() + assert_close(u_hstep.sum(), u) + u += delta_j_B.sum() + assert_close(u_estep.sum(), u) + + +def test_poynting_divergence(sim): + args = {'dxes': sim.dxes, + 'epsilon': sim.epsilon} + + dV = numpy.prod(numpy.meshgrid(*sim.dxes[0], indexing='ij'), axis=0) + + u_eprev = None + for ii in range(1, 8): + u_hstep = fdtd.energy_hstep(e0=sim.es[ii-1], h1=sim.hs[ii], e2=sim.es[ii], **args) + u_estep = fdtd.energy_estep(h0=sim.hs[ii], e1=sim.es[ii], h2=sim.hs[ii + 1], **args) + delta_j_B = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii], dxes=sim.dxes) + + du_half_h2e = u_estep - u_hstep - delta_j_B + div_s_h2e = sim.dt * fdtd.poynting_divergence(e=sim.es[ii], h=sim.hs[ii], dxes=sim.dxes) * dV + assert_fields_close(du_half_h2e, -div_s_h2e, rtol=1e-4) + + if u_eprev is None: + u_eprev = u_estep + continue + + # previous half-step + delta_j_A = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii-1], dxes=sim.dxes) + + du_half_e2h = u_hstep - u_eprev - delta_j_A + div_s_e2h = sim.dt * fdtd.poynting_divergence(e=sim.es[ii-1], h=sim.hs[ii], dxes=sim.dxes) * dV + assert_fields_close(du_half_e2h, -div_s_e2h, rtol=1e-4) + u_eprev = u_estep + + +def test_poynting_planes(sim): + mask = (sim.js[0] != 0) + if mask.sum() > 1: + pytest.skip('test_poynting_planes can only test single point sources') + + args = {'dxes': sim.dxes, + 'epsilon': sim.epsilon} + dV = numpy.prod(numpy.meshgrid(*sim.dxes[0], indexing='ij'), axis=0) + + mx = numpy.roll(mask, (-1, -1), axis=(0, 1)) + my = numpy.roll(mask, -1, axis=2) + mz = numpy.roll(mask, (+1, -1), axis=(0, 3)) + px = numpy.roll(mask, -1, axis=0) + py = mask.copy() + pz = numpy.roll(mask, +1, axis=0) + + u_eprev = None + for ii in range(1, 8): + u_hstep = fdtd.energy_hstep(e0=sim.es[ii-1], h1=sim.hs[ii], e2=sim.es[ii], **args) + u_estep = fdtd.energy_estep(h0=sim.hs[ii], e1=sim.es[ii], h2=sim.hs[ii + 1], **args) + + s_h2e = -fdtd.poynting(e=sim.es[ii], h=sim.hs[ii]) * sim.dt + s_h2e[0] *= sim.dxes[0][1][None, :, None] * sim.dxes[0][2][None, None, :] + s_h2e[1] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][2][None, None, :] + s_h2e[2] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][1][None, :, None] + planes = [s_h2e[px].sum(), -s_h2e[mx].sum(), + s_h2e[py].sum(), -s_h2e[my].sum(), + s_h2e[pz].sum(), -s_h2e[mz].sum()] + assert_close(sum(planes), (u_estep - u_hstep).sum()) + if u_eprev is None: + u_eprev = u_estep + continue + + s_e2h = -fdtd.poynting(e=sim.es[ii - 1], h=sim.hs[ii]) * sim.dt + s_e2h[0] *= sim.dxes[0][1][None, :, None] * sim.dxes[0][2][None, None, :] + s_e2h[1] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][2][None, None, :] + s_e2h[2] *= sim.dxes[0][0][:, None, None] * sim.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()] + assert_close(sum(planes), (u_hstep - u_eprev).sum()) + + # previous half-step + u_eprev = u_estep + +## Now tested elsewhere +#def test_j_dot_e(sim): +# for tt in sim.j_steps: +# e0 = sim.es[tt - 1] +# j1 = sim.js[tt] +# e1 = sim.es[tt] +# +# delta_j_A = fdtd.delta_energy_j(j0=j1, e1=e0, dxes=sim.dxes) +# delta_j_B = fdtd.delta_energy_j(j0=j1, e1=e1, dxes=sim.dxes) +# +# args = {'dxes': sim.dxes, +# 'epsilon': sim.epsilon} +# +# u_eprev = fdtd.energy_estep(h0=sim.hs[tt-1], e1=sim.es[tt-1], h2=sim.hs[tt], **args) +# u_hstep = fdtd.energy_hstep(e0=sim.es[tt-1], h1=sim.hs[tt], e2=sim.es[tt], **args) +# u_estep = fdtd.energy_estep(h0=sim.hs[tt], e1=sim.es[tt], h2=sim.hs[tt + 1], **args) +# +# assert_close(delta_j_A.sum(), (u_hstep - u_eprev).sum(), rtol=1e-4) +# assert_close(delta_j_B.sum(), (u_estep - u_hstep).sum(), rtol=1e-4) + + +@pytest.fixture(scope='module', + params=[(5, 5, 1), + (5, 1, 5), + (5, 5, 5), +# (7, 7, 7), + ]) +def shape(request): + yield (3, *request.param) + + +@pytest.fixture(scope='module', params=[0.3]) +def dt(request): + yield request.param + + +@pytest.fixture(scope='module', params=[1.0, 1.5]) +def epsilon_bg(request): + yield request.param + + +@pytest.fixture(scope='module', params=[1.0, 2.5]) +def epsilon_fg(request): + yield request.param + + +@pytest.fixture(scope='module', params=['center', '000', 'random']) +def epsilon(request, shape, epsilon_bg, epsilon_fg): + is3d = (numpy.array(shape) == 1).sum() == 0 + if is3d: + if request.param == '000': + pytest.skip('Skipping 000 epsilon because test is 3D (for speed)') + if epsilon_bg != 1: + pytest.skip('Skipping epsilon_bg != 1 because test is 3D (for speed)') + if epsilon_fg not in (1.0, 2.0): + pytest.skip('Skipping epsilon_fg not in (1, 2) because test is 3D (for speed)') + + epsilon = numpy.full(shape, epsilon_bg, dtype=float) + if request.param == 'center': + epsilon[:, shape[1]//2, shape[2]//2, shape[3]//2] = epsilon_fg + elif request.param == '000': + epsilon[:, 0, 0, 0] = epsilon_fg + elif request.param == 'random': + epsilon[:] = prng.uniform(low=min(epsilon_bg, epsilon_fg), + high=max(epsilon_bg, epsilon_fg), + size=shape) + + yield epsilon + + +@pytest.fixture(scope='module', params=[1.0])#, 1.5]) +def j_mag(request): + yield request.param + + +@pytest.fixture(scope='module', params=['center', 'random']) +def j_distribution(request, shape, j_mag): + j = numpy.zeros(shape) + if request.param == 'center': + j[:, shape[1]//2, shape[2]//2, shape[3]//2] = j_mag + elif request.param == '000': + j[:, 0, 0, 0] = j_mag + elif request.param == 'random': + j[:] = prng.uniform(low=-j_mag, high=j_mag, size=shape) + yield j + + +@pytest.fixture(scope='module', params=[1.0, 1.5]) +def dx(request): + yield request.param + + +@pytest.fixture(scope='module', params=['uniform']) +def dxes(request, shape, dx): + if request.param == 'uniform': + dxes = [[numpy.full(s, dx) for s in shape[1:]] for _ in range(2)] + yield dxes + + +@pytest.fixture(scope='module', + params=[(0,), + (0, 4, 8), + ] + ) +def j_steps(request): + yield request.param + + +@dataclasses.dataclass() +class SimResult: + shape: Tuple[int] + dt: float + dxes: List[List[numpy.ndarray]] + epsilon: numpy.ndarray + j_distribution: numpy.ndarray + j_steps: Tuple[int] + es: List[numpy.ndarray] = dataclasses.field(default_factory=list) + hs: List[numpy.ndarray] = dataclasses.field(default_factory=list) + js: List[numpy.ndarray] = dataclasses.field(default_factory=list) + + +@pytest.fixture(scope='module') +def sim(request, shape, epsilon, dxes, dt, j_distribution, j_steps): + is3d = (numpy.array(shape) == 1).sum() == 0 + if is3d: + if dt != 0.3: + pytest.skip('Skipping dt != 0.3 because test is 3D (for speed)') + + sim = SimResult( + shape=shape, + dt=dt, + dxes=dxes, + epsilon=epsilon, + j_distribution=j_distribution, + j_steps=j_steps, + ) + + e = numpy.zeros_like(epsilon) + h = numpy.zeros_like(epsilon) + + assert 0 in j_steps + j_zeros = numpy.zeros_like(j_distribution) + + eh2h = fdtd.maxwell_h(dt=dt, dxes=dxes) + eh2e = fdtd.maxwell_e(dt=dt, dxes=dxes) + for tt in range(10): + e = e.copy() + h = h.copy() + eh2h(e, h) + eh2e(e, h, epsilon) + if tt in j_steps: + e += j_distribution / epsilon + sim.js.append(j_distribution) + else: + sim.js.append(j_zeros) + sim.es.append(e) + sim.hs.append(h) + return sim + + diff --git a/fdfd_tools/test_fdtd.py b/fdfd_tools/test_fdtd.py deleted file mode 100644 index 15c5b4d..0000000 --- a/fdfd_tools/test_fdtd.py +++ /dev/null @@ -1,353 +0,0 @@ -import unittest -import numpy - -from fdfd_tools import fdtd - - -class BasicTests(): - def test_initial_fields(self): - # Make sure initial fields didn't change - e0 = self.es[0] - h0 = self.hs[0] - mask = self.src_mask - - self.assertEqual(e0[mask], self.j_mag / self.epsilon[mask]) - self.assertFalse(e0[~mask].any()) - self.assertFalse(h0.any()) - - - def test_initial_energy(self): - e0 = self.es[0] - h0 = self.hs[0] - h1 = self.hs[1] - mask = self.src_mask[1] - dxes = self.dxes if self.dxes is not None else tuple(tuple(numpy.ones(s) for s in e0.shape[1:]) for _ in range(2)) - dV = numpy.prod(numpy.meshgrid(*dxes[0], indexing='ij'), axis=0) - u0 = self.j_mag * self.j_mag / self.epsilon[self.src_mask] * dV[mask] - args = {'dxes': self.dxes, - 'epsilon': self.epsilon} - - # Make sure initial energy and E dot J are correct - energy0 = fdtd.energy_estep(h0=h0, e1=e0, h2=self.hs[1], **args) - e_dot_j_0 = fdtd.delta_energy_j(j0=(e0 - 0) * self.epsilon, e1=e0, dxes=self.dxes) - self.assertEqual(energy0[mask], u0) - self.assertFalse(energy0[~mask].any(), msg='energy0: {}'.format(energy0)) - self.assertEqual(e_dot_j_0[mask], u0) - self.assertFalse(e_dot_j_0[~mask].any(), msg='e_dot_j_0: {}'.format(e_dot_j_0)) - - - def test_energy_conservation(self): - e0 = self.es[0] - u0 = fdtd.delta_energy_j(j0=(e0 - 0) * self.epsilon, e1=e0, dxes=self.dxes).sum() - args = {'dxes': self.dxes, - 'epsilon': self.epsilon} - - 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) - self.assertTrue(numpy.allclose(u_hstep.sum(), u0), msg='u_hstep: {}\n{}'.format(u_hstep.sum(), numpy.rollaxis(u_hstep, -1))) - self.assertTrue(numpy.allclose(u_estep.sum(), u0), msg='u_estep: {}\n{}'.format(u_estep.sum(), numpy.rollaxis(u_estep, -1))) - - - def test_poynting_divergence(self): - args = {'dxes': self.dxes, - 'epsilon': self.epsilon} - - 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) - - 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) - - du_half_h2e = u_estep - u_hstep - div_s_h2e = self.dt * fdtd.poynting_divergence(e=self.es[ii], h=self.hs[ii], dxes=self.dxes) * dV - self.assertTrue(numpy.allclose(du_half_h2e, -div_s_h2e, rtol=1e-4), - msg='du_half_h2e\n{}\ndiv_s_h2e\n{}'.format(numpy.rollaxis(du_half_h2e, -1), - -numpy.rollaxis(div_s_h2e, -1))) - - if u_eprev is None: - u_eprev = u_estep - continue - - # previous half-step - du_half_e2h = u_hstep - u_eprev - div_s_e2h = self.dt * fdtd.poynting_divergence(e=self.es[ii-1], h=self.hs[ii], dxes=self.dxes) * dV - self.assertTrue(numpy.allclose(du_half_e2h, -div_s_e2h, rtol=1e-4), - msg='du_half_e2h\n{}\ndiv_s_e2h\n{}'.format(numpy.rollaxis(du_half_e2h, -1), - -numpy.rollaxis(div_s_e2h, -1))) - u_eprev = u_estep - - - def test_poynting_planes(self): - args = {'dxes': self.dxes, - 'epsilon': self.epsilon} - 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) - - 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, :] - s_h2e[2] *= dxes[0][0][:, None, None] * dxes[0][1][None, :, None] - planes = [s_h2e[px].sum(), -s_h2e[mx].sum(), - s_h2e[py].sum(), -s_h2e[my].sum(), - s_h2e[pz].sum(), -s_h2e[mz].sum()] - 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): - shape = [3, 5, 5, 1] - self.dt = 0.5 - self.epsilon = numpy.ones(shape, dtype=float) - self.j_mag = 32 - self.dxes = None - - self.src_mask = numpy.zeros_like(self.epsilon, dtype=bool) - self.src_mask[1, 2, 2, 0] = True - - e = numpy.zeros_like(self.epsilon) - h = numpy.zeros_like(self.epsilon) - e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask] - self.es = [e] - self.hs = [h] - - eh2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes) - eh2e = fdtd.maxwell_e(dt=self.dt, dxes=self.dxes) - for _ in range(9): - e = e.copy() - h = h.copy() - eh2h(e, h) - eh2e(e, h, self.epsilon) - self.es.append(e) - self.hs.append(h) - - -class Basic2DUniformDX3(unittest.TestCase, BasicTests): - def setUp(self): - shape = [3, 5, 5, 1] - self.dt = 0.5 - self.j_mag = 32 - self.dxes = tuple(tuple(numpy.full(s, 2.0) for s in shape[1:]) for _ in range(2)) - - self.src_mask = numpy.zeros(shape, dtype=bool) - self.src_mask[1, 2, 2, 0] = True - - self.epsilon = numpy.full(shape, 1, dtype=float) - self.epsilon[self.src_mask] = 2 - - e = numpy.zeros_like(self.epsilon) - h = numpy.zeros_like(self.epsilon) - e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask] - self.es = [e] - self.hs = [h] - - eh2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes) - eh2e = fdtd.maxwell_e(dt=self.dt, dxes=self.dxes) - for _ in range(9): - e = e.copy() - h = h.copy() - eh2h(e, h) - eh2e(e, h, self.epsilon) - self.es.append(e) - self.hs.append(h) - - -class Basic3DUniformDXOnlyVacuum(unittest.TestCase, BasicTests): - def setUp(self): - shape = [3, 5, 5, 5] - self.dt = 0.5 - self.epsilon = numpy.ones(shape, dtype=float) - self.j_mag = 32 - self.dxes = tuple(tuple(numpy.ones(s) for s in shape[1:]) for _ in range(2)) - - self.src_mask = numpy.zeros_like(self.epsilon, dtype=bool) - self.src_mask[1, 2, 2, 2] = True - - e = numpy.zeros_like(self.epsilon) - h = numpy.zeros_like(self.epsilon) - e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask] - self.es = [e] - self.hs = [h] - - eh2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes) - eh2e = fdtd.maxwell_e(dt=self.dt, dxes=self.dxes) - for _ in range(9): - e = e.copy() - h = h.copy() - eh2h(e, h) - eh2e(e, h, self.epsilon) - self.es.append(e) - self.hs.append(h) - - - -class Basic3DUniformDXUniformN(unittest.TestCase, BasicTests): - def setUp(self): - shape = [3, 5, 5, 5] - self.dt = 0.5 - self.epsilon = numpy.full(shape, 2, dtype=float) - self.j_mag = 32 - self.dxes = tuple(tuple(numpy.ones(s) for s in shape[1:]) for _ in range(2)) - - self.src_mask = numpy.zeros_like(self.epsilon, dtype=bool) - self.src_mask[1, 2, 2, 2] = True - - e = numpy.zeros_like(self.epsilon) - h = numpy.zeros_like(self.epsilon) - e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask] - self.es = [e] - self.hs = [h] - - eh2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes) - eh2e = fdtd.maxwell_e(dt=self.dt, dxes=self.dxes) - for _ in range(9): - e = e.copy() - h = h.copy() - eh2h(e, h) - eh2e(e, h, self.epsilon) - self.es.append(e) - self.hs.append(h) - - -class Basic3DUniformDX(unittest.TestCase, BasicTests): - def setUp(self): - shape = [3, 5, 5, 5] - self.dt = 0.33 - self.j_mag = 32 - self.dxes = tuple(tuple(numpy.ones(s) for s in shape[1:]) for _ in range(2)) - - self.src_mask = numpy.zeros(shape, dtype=bool) - self.src_mask[1, 2, 2, 2] = True - - self.epsilon = numpy.full(shape, 1, dtype=float) - self.epsilon[self.src_mask] = 2 - - e = numpy.zeros_like(self.epsilon) - h = numpy.zeros_like(self.epsilon) - e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask] - self.es = [e] - self.hs = [h] - - eh2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes) - eh2e = fdtd.maxwell_e(dt=self.dt, dxes=self.dxes) - for _ in range(9): - e = e.copy() - h = h.copy() - eh2h(e, h) - eh2e(e, h, self.epsilon) - self.es.append(e) - self.hs.append(h) - - -class Basic3DUniformDX3(unittest.TestCase, BasicTests): - def setUp(self): - shape = [3, 5, 5, 5] - self.dt = 0.5 - self.j_mag = 32 - self.dxes = tuple(tuple(numpy.full(s, 3.0) for s in shape[1:]) for _ in range(2)) - - self.src_mask = numpy.zeros(shape, dtype=bool) - self.src_mask[1, 2, 2, 2] = True - - self.epsilon = numpy.full(shape, 1, dtype=float) - self.epsilon[self.src_mask] = 2 - - e = numpy.zeros_like(self.epsilon) - h = numpy.zeros_like(self.epsilon) - e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask] - self.es = [e] - self.hs = [h] - - eh2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes) - eh2e = fdtd.maxwell_e(dt=self.dt, dxes=self.dxes) - for _ in range(9): - e = e.copy() - h = h.copy() - eh2h(e, h) - eh2e(e, h, self.epsilon) - self.es.append(e) - self.hs.append(h) - - -class JdotE_3DUniformDX(unittest.TestCase): - def setUp(self): - shape = [3, 5, 5, 5] - self.dt = 0.5 - self.j_mag = 32 - self.dxes = tuple(tuple(numpy.full(s, 2.0) for s in shape[1:]) for _ in range(2)) - - self.src_mask = numpy.zeros(shape, dtype=bool) - self.src_mask[1, 2, 2, 2] = True - - self.epsilon = numpy.full(shape, 4, dtype=float) - self.epsilon[self.src_mask] = 2 - - e = numpy.random.randint(-128, 128 + 1, size=shape).astype(float) - h = numpy.random.randint(-128, 128 + 1, size=shape).astype(float) - self.es = [e] - self.hs = [h] - - eh2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes) - eh2e = fdtd.maxwell_e(dt=self.dt, dxes=self.dxes) - for ii in range(9): - e = e.copy() - h = h.copy() - eh2h(e, h) - eh2e(e, h, self.epsilon) - self.es.append(e) - self.hs.append(h) - - if ii == 1: - e[self.src_mask] += self.j_mag / self.epsilon[self.src_mask] - self.j_dot_e = self.j_mag * e[self.src_mask] - - - def test_j_dot_e(self): - e0 = self.es[2] - j0 = numpy.zeros_like(e0) - j0[self.src_mask] = self.j_mag - u0 = fdtd.delta_energy_j(j0=j0, e1=e0, dxes=self.dxes) - args = {'dxes': self.dxes, - 'epsilon': self.epsilon} - - ii=2 - 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) - #print(u0.sum(), (u_estep - u_hstep).sum()) - self.assertTrue(numpy.allclose(u0.sum(), (u_estep - u_hstep).sum(), rtol=1e-4)) -