From 2a1fa1df7bdf79c285fe82d655bd58b01d24c3d4 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 27 Oct 2019 12:43:06 -0700 Subject: [PATCH] Take care of dxes at the poynting() level rather than poynting_divergence --- meanas/fdtd/energy.py | 34 +++++++++++++++++++++------------- meanas/test/test_fdtd.py | 18 +++++------------- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/meanas/fdtd/energy.py b/meanas/fdtd/energy.py index 26fb036..ffbbf80 100644 --- a/meanas/fdtd/energy.py +++ b/meanas/fdtd/energy.py @@ -4,23 +4,31 @@ import numpy from .. import dx_lists_t, field_t, field_updater -def poynting(e, h): - s = (numpy.roll(e[1], -1, axis=0) * h[2] - numpy.roll(e[2], -1, axis=0) * h[1], - numpy.roll(e[2], -1, axis=1) * h[0] - numpy.roll(e[0], -1, axis=1) * h[2], - numpy.roll(e[0], -1, axis=2) * h[1] - numpy.roll(e[1], -1, axis=2) * h[0]) - return numpy.array(s) - - -def poynting_divergence(s=None, *, e=None, h=None, dxes=None): # TODO dxes +def poynting(e, h, dxes=None): if dxes is None: dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2)) - if s is None: - s = poynting(e, h) + ex = e[0] * dxes[0][0][:, None, None] + ey = e[1] * dxes[0][1][None, :, None] + ez = e[2] * dxes[0][2][None, None, :] + hx = h[0] * dxes[1][0][:, None, None] + hy = h[1] * dxes[1][1][None, :, None] + hz = h[2] * dxes[1][2][None, None, :] - ds = ((s[0] - numpy.roll(s[0], 1, axis=0)) / numpy.sqrt(dxes[0][0] * dxes[1][0])[:, None, None] + - (s[1] - numpy.roll(s[1], 1, axis=1)) / numpy.sqrt(dxes[0][1] * dxes[1][1])[None, :, None] + - (s[2] - numpy.roll(s[2], 1, axis=2)) / numpy.sqrt(dxes[0][2] * dxes[1][2])[None, None, :] ) + s = numpy.empty_like(e) + s[0] = numpy.roll(ey, -1, axis=0) * hz - numpy.roll(ez, -1, axis=0) * hy + s[1] = numpy.roll(ez, -1, axis=1) * hx - numpy.roll(ex, -1, axis=1) * hz + s[2] = numpy.roll(ex, -1, axis=2) * hy - numpy.roll(ey, -1, axis=2) * hx + return s + + +def poynting_divergence(s=None, *, e=None, h=None, dxes=None): # TODO dxes + if s is None: + s = poynting(e, h, dxes=dxes) + + ds = ((s[0] - numpy.roll(s[0], 1, axis=0)) + + (s[1] - numpy.roll(s[1], 1, axis=1)) + + (s[2] - numpy.roll(s[2], 1, axis=2)) ) return ds diff --git a/meanas/test/test_fdtd.py b/meanas/test/test_fdtd.py index 08b678e..1b185f5 100644 --- a/meanas/test/test_fdtd.py +++ b/meanas/test/test_fdtd.py @@ -77,16 +77,14 @@ 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) + 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 + div_s_h2e = sim.dt * fdtd.poynting_divergence(e=sim.es[ii], h=sim.hs[ii], dxes=sim.dxes) assert_fields_close(du_half_h2e, -div_s_h2e) if u_eprev is None: @@ -97,7 +95,7 @@ def test_poynting_divergence(sim): 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 + div_s_e2h = sim.dt * fdtd.poynting_divergence(e=sim.es[ii-1], h=sim.hs[ii], dxes=sim.dxes) assert_fields_close(du_half_e2h, -div_s_e2h) u_eprev = u_estep @@ -123,10 +121,7 @@ def test_poynting_planes(sim): 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] + s_h2e = -fdtd.poynting(e=sim.es[ii], h=sim.hs[ii], dxes=sim.dxes) * sim.dt 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()] @@ -135,10 +130,7 @@ def test_poynting_planes(sim): 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] + s_e2h = -fdtd.poynting(e=sim.es[ii - 1], h=sim.hs[ii], dxes=sim.dxes) * sim.dt 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()]