jan/meanas
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Take care of dxes at the poynting() level rather than poynting_divergence

Browse Source
This commit is contained in:
Jan Petykiewicz 2019-10-27 12:43:06 -07:00
parent c0f14cd07b
commit 2a1fa1df7b
2 changed files with 26 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
meanas/fdtd/energy.py
View File

@ -4,23 +4,31 @@ import numpy
from .. import dx_lists_t, field_t, field_updater from .. import dx_lists_t, field_t, field_updater
def poynting(e, h): def poynting(e, h, dxes=None):
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
if dxes is None: if dxes is None:
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2)) dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
if s is None: ex = e[0] * dxes[0][0][:, None, None]
s = poynting(e, h) 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 = numpy.empty_like(e)
(s[1] - numpy.roll(s[1], 1, axis=1)) / numpy.sqrt(dxes[0][1] * dxes[1][1])[None, :, None] + s[0] = numpy.roll(ey, -1, axis=0) * hz - numpy.roll(ez, -1, axis=0) * hy
(s[2] - numpy.roll(s[2], 1, axis=2)) / numpy.sqrt(dxes[0][2] * dxes[1][2])[None, None, :] ) 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 return ds

18
meanas/test/test_fdtd.py
View File

@ -77,16 +77,14 @@ def test_poynting_divergence(sim):
args = {'dxes': sim.dxes, args = {'dxes': sim.dxes,
'epsilon': sim.epsilon} 'epsilon': sim.epsilon}
dV = numpy.prod(numpy.meshgrid(*sim.dxes[0], indexing='ij'), axis=0)
u_eprev = None u_eprev = None
for ii in range(1, 8): 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_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) 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 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) assert_fields_close(du_half_h2e, -div_s_h2e)
if u_eprev is None: 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) 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 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) assert_fields_close(du_half_e2h, -div_s_e2h)
u_eprev = u_estep 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_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) 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 = -fdtd.poynting(e=sim.es[ii], h=sim.hs[ii], dxes=sim.dxes) * 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(), planes = [s_h2e[px].sum(), -s_h2e[mx].sum(),
s_h2e[py].sum(), -s_h2e[my].sum(), s_h2e[py].sum(), -s_h2e[my].sum(),
s_h2e[pz].sum(), -s_h2e[mz].sum()] s_h2e[pz].sum(), -s_h2e[mz].sum()]
@ -135,10 +130,7 @@ def test_poynting_planes(sim):
u_eprev = u_estep u_eprev = u_estep
continue continue
s_e2h = -fdtd.poynting(e=sim.es[ii - 1], h=sim.hs[ii]) * sim.dt s_e2h = -fdtd.poynting(e=sim.es[ii - 1], h=sim.hs[ii], dxes=sim.dxes) * 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(), planes = [s_e2h[px].sum(), -s_e2h[mx].sum(),
s_e2h[py].sum(), -s_e2h[my].sum(), s_e2h[py].sum(), -s_e2h[my].sum(),
s_e2h[pz].sum(), -s_e2h[mz].sum()] s_e2h[pz].sum(), -s_e2h[mz].sum()]