333 lines
11 KiB
Python
333 lines
11 KiB
Python
import numpy
|
|
|
|
from ..fdmath import dx_lists_t, fdfield_t, fdfield
|
|
from ..fdmath.functional import deriv_back
|
|
|
|
|
|
"""
|
|
Energy- and flux-accounting helpers for Yee-grid FDTD fields.
|
|
|
|
These functions complement the derivation in `meanas.fdtd`:
|
|
|
|
- `poynting(...)` and `poynting_divergence(...)` evaluate the discrete flux terms
|
|
from the exact time-domain Poynting identity.
|
|
- `energy_hstep(...)` / `energy_estep(...)` evaluate the two staggered energy
|
|
expressions.
|
|
- `delta_energy_*` helpers evaluate the corresponding energy changes between
|
|
adjacent half-steps.
|
|
|
|
The helpers are intended for diagnostics, regression tests, and consistency
|
|
checks between source work, field energy, and flux through cell faces.
|
|
"""
|
|
|
|
|
|
def poynting(
|
|
e: fdfield,
|
|
h: fdfield,
|
|
dxes: dx_lists_t | None = None,
|
|
) -> fdfield_t:
|
|
r"""
|
|
Calculate the poynting vector `S` ($S$).
|
|
|
|
This is the energy transfer rate (amount of energy `U` per `dt` transferred
|
|
between adjacent cells) in each direction that happens during the half-step
|
|
bounded by the two provided fields.
|
|
|
|
The returned vector field `S` is the energy flow across +x, +y, and +z
|
|
boundaries of the corresponding `U` cell. For example,
|
|
|
|
```
|
|
mx = numpy.roll(mask, -1, axis=0)
|
|
my = numpy.roll(mask, -1, axis=1)
|
|
mz = numpy.roll(mask, -1, axis=2)
|
|
|
|
u_hstep = fdtd.energy_hstep(e0=es[ii - 1], h1=hs[ii], e2=es[ii], **args)
|
|
u_estep = fdtd.energy_estep(h0=hs[ii], e1=es[ii], h2=hs[ii + 1], **args)
|
|
delta_j_B = fdtd.delta_energy_j(j0=js[ii], e1=es[ii], dxes=dxes)
|
|
du_half_h2e = u_estep - u_hstep - delta_j_B
|
|
|
|
s_h2e = -fdtd.poynting(e=es[ii], h=hs[ii], dxes=dxes) * dt
|
|
planes = [s_h2e[0, mask].sum(), -s_h2e[0, mx].sum(),
|
|
s_h2e[1, mask].sum(), -s_h2e[1, my].sum(),
|
|
s_h2e[2, mask].sum(), -s_h2e[2, mz].sum()]
|
|
|
|
assert_close(sum(planes), du_half_h2e[mask])
|
|
```
|
|
|
|
(see `meanas.tests.test_fdtd.test_poynting_planes`)
|
|
|
|
The full relationship is
|
|
$$
|
|
\begin{aligned}
|
|
(U_{l+\frac{1}{2}} - U_l) / \Delta_t
|
|
&= -\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} \\
|
|
- \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\
|
|
- \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
|
|
(U_l - U_{l-\frac{1}{2}}) / \Delta_t
|
|
&= -\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} \\
|
|
- \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\
|
|
- \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
|
|
\end{aligned}
|
|
$$
|
|
|
|
These equalities are exact and should practically hold to within numerical precision.
|
|
No time- or spatial-averaging is necessary. (See `meanas.fdtd` docs for derivation.)
|
|
|
|
Args:
|
|
e: E-field
|
|
h: H-field (one half-timestep before or after `e`)
|
|
dxes: Grid description; see `meanas.fdmath`.
|
|
|
|
Returns:
|
|
s: Vector field. Components indicate the energy transfer rate from the
|
|
corresponding energy cell into its +x, +y, and +z neighbors during
|
|
the half-step from the time of the earlier input field until the
|
|
time of later input field.
|
|
"""
|
|
if dxes is None:
|
|
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
|
|
|
|
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, :]
|
|
|
|
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 fdfield_t(s)
|
|
|
|
|
|
def poynting_divergence(
|
|
s: fdfield | None = None,
|
|
*,
|
|
e: fdfield | None = None,
|
|
h: fdfield | None = None,
|
|
dxes: dx_lists_t | None = None,
|
|
) -> fdfield_t:
|
|
"""
|
|
Calculate the divergence of the poynting vector.
|
|
|
|
This is the net energy flow for each cell, i.e. the change in energy `U`
|
|
per `dt` caused by transfer of energy to nearby cells (rather than
|
|
absorption/emission by currents `J` or `M`).
|
|
|
|
See `poynting` and `meanas.fdtd` for more details.
|
|
Args:
|
|
s: Poynting vector, as calculated with `poynting`. Optional; caller
|
|
can provide `e` and `h` instead.
|
|
e: E-field (optional; need either `s` or both `e` and `h`)
|
|
h: H-field (optional; need either `s` or both `e` and `h`)
|
|
dxes: Grid description; see `meanas.fdmath`.
|
|
|
|
Returns:
|
|
ds: Divergence of the poynting vector.
|
|
Entries indicate the net energy flow out of the corresponding
|
|
energy cell.
|
|
"""
|
|
if s is None:
|
|
assert e is not None
|
|
assert h is not None
|
|
assert dxes is not None
|
|
s = poynting(e, h, dxes=dxes)
|
|
|
|
Dx, Dy, Dz = deriv_back()
|
|
ds = Dx(s[0]) + Dy(s[1]) + Dz(s[2])
|
|
return fdfield_t(ds)
|
|
|
|
|
|
def energy_hstep(
|
|
e0: fdfield,
|
|
h1: fdfield,
|
|
e2: fdfield,
|
|
epsilon: fdfield | None = None,
|
|
mu: fdfield | None = None,
|
|
dxes: dx_lists_t | None = None,
|
|
) -> fdfield_t:
|
|
"""
|
|
Calculate energy `U` at the time of the provided H-field `h1`.
|
|
|
|
TODO: Figure out what this means spatially.
|
|
|
|
Args:
|
|
e0: E-field one half-timestep before the energy.
|
|
h1: H-field (at the same timestep as the energy).
|
|
e2: E-field one half-timestep after the energy.
|
|
epsilon: Dielectric constant distribution.
|
|
mu: Magnetic permeability distribution.
|
|
dxes: Grid description; see `meanas.fdmath`.
|
|
|
|
Returns:
|
|
Energy, at the time of the H-field `h1`.
|
|
"""
|
|
u = dxmul(e0 * e2, h1 * h1, epsilon, mu, dxes)
|
|
return fdfield_t(u)
|
|
|
|
|
|
def energy_estep(
|
|
h0: fdfield,
|
|
e1: fdfield,
|
|
h2: fdfield,
|
|
epsilon: fdfield | None = None,
|
|
mu: fdfield | None = None,
|
|
dxes: dx_lists_t | None = None,
|
|
) -> fdfield_t:
|
|
"""
|
|
Calculate energy `U` at the time of the provided E-field `e1`.
|
|
|
|
TODO: Figure out what this means spatially.
|
|
|
|
Args:
|
|
h0: H-field one half-timestep before the energy.
|
|
e1: E-field (at the same timestep as the energy).
|
|
h2: H-field one half-timestep after the energy.
|
|
epsilon: Dielectric constant distribution.
|
|
mu: Magnetic permeability distribution.
|
|
dxes: Grid description; see `meanas.fdmath`.
|
|
|
|
Returns:
|
|
Energy, at the time of the E-field `e1`.
|
|
"""
|
|
u = dxmul(e1 * e1, h0 * h2, epsilon, mu, dxes)
|
|
return fdfield_t(u)
|
|
|
|
|
|
def delta_energy_h2e(
|
|
dt: float,
|
|
e0: fdfield,
|
|
h1: fdfield,
|
|
e2: fdfield,
|
|
h3: fdfield,
|
|
epsilon: fdfield | None = None,
|
|
mu: fdfield | None = None,
|
|
dxes: dx_lists_t | None = None,
|
|
) -> fdfield_t:
|
|
"""
|
|
Change in energy during the half-step from `h1` to `e2`.
|
|
|
|
This is just from (e2 * e2 + h3 * h1) - (h1 * h1 + e0 * e2)
|
|
|
|
Args:
|
|
e0: E-field one half-timestep before the start of the energy delta.
|
|
h1: H-field at the start of the energy delta.
|
|
e2: E-field at the end of the energy delta (one half-timestep after `h1`).
|
|
h3: H-field one half-timestep after the end of the energy delta.
|
|
epsilon: Dielectric constant distribution.
|
|
mu: Magnetic permeability distribution.
|
|
dxes: Grid description; see `meanas.fdmath`.
|
|
|
|
Returns:
|
|
Change in energy from the time of `h1` to the time of `e2`.
|
|
"""
|
|
de = e2 * (e2 - e0) / dt
|
|
dh = h1 * (h3 - h1) / dt
|
|
du = dxmul(de, dh, epsilon, mu, dxes)
|
|
return fdfield_t(du)
|
|
|
|
|
|
def delta_energy_e2h(
|
|
dt: float,
|
|
h0: fdfield,
|
|
e1: fdfield,
|
|
h2: fdfield,
|
|
e3: fdfield,
|
|
epsilon: fdfield | None = None,
|
|
mu: fdfield | None = None,
|
|
dxes: dx_lists_t | None = None,
|
|
) -> fdfield_t:
|
|
"""
|
|
Change in energy during the half-step from `e1` to `h2`.
|
|
|
|
This is just from (h2 * h2 + e3 * e1) - (e1 * e1 + h0 * h2)
|
|
|
|
Args:
|
|
h0: E-field one half-timestep before the start of the energy delta.
|
|
e1: H-field at the start of the energy delta.
|
|
h2: E-field at the end of the energy delta (one half-timestep after `e1`).
|
|
e3: H-field one half-timestep after the end of the energy delta.
|
|
epsilon: Dielectric constant distribution.
|
|
mu: Magnetic permeability distribution.
|
|
dxes: Grid description; see `meanas.fdmath`.
|
|
|
|
Returns:
|
|
Change in energy from the time of `e1` to the time of `h2`.
|
|
"""
|
|
de = e1 * (e3 - e1) / dt
|
|
dh = h2 * (h2 - h0) / dt
|
|
du = dxmul(de, dh, epsilon, mu, dxes)
|
|
return fdfield_t(du)
|
|
|
|
|
|
def delta_energy_j(
|
|
j0: fdfield,
|
|
e1: fdfield,
|
|
dxes: dx_lists_t | None = None,
|
|
) -> fdfield_t:
|
|
r"""
|
|
Calculate the electric-current work term $J \cdot E$ on the Yee grid.
|
|
|
|
This is the source contribution that appears beside the flux divergence in
|
|
the discrete Poynting identities documented in `meanas.fdtd`.
|
|
|
|
Note that each value of `J` contributes twice in a full Yee cycle (once per
|
|
half-step energy balance) even though it directly changes `E` only once.
|
|
|
|
Args:
|
|
j0: Electric-current density sampled at the same half-step as the
|
|
current work term.
|
|
e1: Electric field sampled at the matching integer timestep.
|
|
dxes: Grid description; defaults to unit spacing.
|
|
|
|
Returns:
|
|
Per-cell source-work contribution with the scalar field shape.
|
|
"""
|
|
if dxes is None:
|
|
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
|
|
|
|
du = ((j0 * e1).sum(axis=0)
|
|
* dxes[0][0][:, None, None]
|
|
* dxes[0][1][None, :, None]
|
|
* dxes[0][2][None, None, :])
|
|
return fdfield_t(du)
|
|
|
|
|
|
def dxmul(
|
|
ee: fdfield,
|
|
hh: fdfield,
|
|
epsilon: fdfield | float | None = None,
|
|
mu: fdfield | float | None = None,
|
|
dxes: dx_lists_t | None = None,
|
|
) -> fdfield_t:
|
|
"""
|
|
Multiply E- and H-like field products by material weights and cell volumes.
|
|
|
|
Args:
|
|
ee: Three-component electric-field product, such as `e0 * e2`.
|
|
hh: Three-component magnetic-field product, such as `h1 * h1`.
|
|
epsilon: Electric material weight; defaults to `1`.
|
|
mu: Magnetic material weight; defaults to `1`.
|
|
dxes: Grid description; defaults to unit spacing.
|
|
|
|
Returns:
|
|
Scalar field containing the weighted electric plus magnetic contribution
|
|
for each Yee cell.
|
|
"""
|
|
if epsilon is None:
|
|
epsilon = 1
|
|
if mu is None:
|
|
mu = 1
|
|
if dxes is None:
|
|
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
|
|
|
|
result = ((ee * epsilon).sum(axis=0)
|
|
* dxes[0][0][:, None, None]
|
|
* dxes[0][1][None, :, None]
|
|
* dxes[0][2][None, None, :]
|
|
+ (hh * mu).sum(axis=0)
|
|
* dxes[1][0][:, None, None]
|
|
* dxes[1][1][None, :, None]
|
|
* dxes[1][2][None, None, :])
|
|
return fdfield_t(result)
|