274 lines
11 KiB
Python
274 lines
11 KiB
Python
r"""
|
|
Utilities for running finite-difference time-domain (FDTD) simulations
|
|
|
|
See the discussion of `Maxwell's Equations` in `meanas.fdmath` for basic
|
|
mathematical background.
|
|
|
|
|
|
Timestep
|
|
========
|
|
|
|
From the discussion of "Plane waves and the Dispersion relation" in `meanas.fdmath`,
|
|
we have
|
|
|
|
$$ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} < 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) $$
|
|
|
|
or, if $\Delta_x = \Delta_y = \Delta_z$, then $c \Delta_t < \frac{\Delta_x}{\sqrt{3}}$.
|
|
|
|
Based on this, we can set
|
|
|
|
dt = sqrt(mu.min() * epsilon.min()) / sqrt(1/dx_min**2 + 1/dy_min**2 + 1/dz_min**2)
|
|
|
|
The `dx_min`, `dy_min`, `dz_min` should be the minimum value across both the base and derived grids.
|
|
|
|
|
|
Poynting Vector and Energy Conservation
|
|
=======================================
|
|
|
|
Let
|
|
|
|
$$ \begin{aligned}
|
|
\tilde{S}_{l, l', \vec{r}} &=& &\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}} \\
|
|
&=& &\vec{x} (\tilde{E}^y_{l,m+1,n,p} \hat{H}^z_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^z_{l,m+1,n,p} \hat{H}^y_{l', \vec{r} + \frac{1}{2}}) \\
|
|
& &+ &\vec{y} (\tilde{E}^z_{l,m,n+1,p} \hat{H}^x_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^x_{l,m,n+1,p} \hat{H}^z_{l', \vec{r} + \frac{1}{2}}) \\
|
|
& &+ &\vec{z} (\tilde{E}^x_{l,m,n,p+1} \hat{H}^y_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^y_{l,m,n,p+1} \hat{H}^z_{l', \vec{r} + \frac{1}{2}})
|
|
\end{aligned}
|
|
$$
|
|
|
|
where $\vec{r} = (m, n, p)$ and $\otimes$ is a modified cross product
|
|
in which the $\tilde{E}$ terms are shifted as indicated.
|
|
|
|
By taking the divergence and rearranging terms, we can show that
|
|
|
|
$$
|
|
\begin{aligned}
|
|
\hat{\nabla} \cdot \tilde{S}_{l, l', \vec{r}}
|
|
&= \hat{\nabla} \cdot (\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}}) \\
|
|
&= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l, \vec{r}} -
|
|
\tilde{E}_{l, \vec{r}} \cdot \hat{\nabla} \times \hat{H}_{l', \vec{r} + \frac{1}{2}} \\
|
|
&= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot
|
|
(-\tilde{\partial}_t \mu_{\vec{r} + \frac{1}{2}} \hat{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} -
|
|
\hat{M}_{l, \vec{r} + \frac{1}{2}}) -
|
|
\tilde{E}_{l, \vec{r}} \cdot (\hat{\partial}_t \tilde{\epsilon}_{\vec{r}} \tilde{E}_{l'+\frac{1}{2}, \vec{r}} +
|
|
\tilde{J}_{l', \vec{r}}) \\
|
|
&= \hat{H}_{l'} \cdot (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) -
|
|
\tilde{E}_l \cdot (\epsilon / \Delta_t )(\tilde{E}_{l'+\frac{1}{2}} - \tilde{E}_{l'-\frac{1}{2}})
|
|
- \hat{H}_{l'} \cdot \hat{M}_{l} - \tilde{E}_l \cdot \tilde{J}_{l'} \\
|
|
\end{aligned}
|
|
$$
|
|
|
|
where in the last line the spatial subscripts have been dropped to emphasize
|
|
the time subscripts $l, l'$, i.e.
|
|
|
|
$$
|
|
\begin{aligned}
|
|
\tilde{E}_l &= \tilde{E}_{l, \vec{r}} \\
|
|
\hat{H}_l &= \tilde{H}_{l, \vec{r} + \frac{1}{2}} \\
|
|
\tilde{\epsilon} &= \tilde{\epsilon}_{\vec{r}} \\
|
|
\end{aligned}
|
|
$$
|
|
|
|
etc.
|
|
For $l' = l + \frac{1}{2}$ we get
|
|
|
|
$$
|
|
\begin{aligned}
|
|
\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}}
|
|
&= \hat{H}_{l + \frac{1}{2}} \cdot
|
|
(-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) -
|
|
\tilde{E}_l \cdot (\epsilon / \Delta_t)(\tilde{E}_{l+1} - \tilde{E}_l)
|
|
- \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\
|
|
&= (-\mu / \Delta_t)(\hat{H}^2_{l + \frac{1}{2}} - \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}) -
|
|
(\epsilon / \Delta_t)(\tilde{E}_{l+1} \cdot \tilde{E}_l - \tilde{E}^2_l)
|
|
- \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\
|
|
&= -(\mu \hat{H}^2_{l + \frac{1}{2}}
|
|
+\epsilon \tilde{E}_{l+1} \cdot \tilde{E}_l) / \Delta_t \\
|
|
+(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}
|
|
+\epsilon \tilde{E}^2_l) / \Delta_t \\
|
|
- \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\
|
|
- \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
|
|
\end{aligned}
|
|
$$
|
|
|
|
and for $l' = l - \frac{1}{2}$,
|
|
|
|
$$
|
|
\begin{aligned}
|
|
\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}}
|
|
&= (\mu \hat{H}^2_{l - \frac{1}{2}}
|
|
+\epsilon \tilde{E}_{l-1} \cdot \tilde{E}_l) / \Delta_t \\
|
|
-(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}
|
|
+\epsilon \tilde{E}^2_l) / \Delta_t \\
|
|
- \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\
|
|
- \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
|
|
\end{aligned}
|
|
$$
|
|
|
|
These two results form the discrete time-domain analogue to Poynting's theorem.
|
|
They hint at the expressions for the energy, which can be calculated at the same
|
|
time-index as either the E or H field:
|
|
|
|
$$
|
|
\begin{aligned}
|
|
U_l &= \epsilon \tilde{E}^2_l + \mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}} \\
|
|
U_{l + \frac{1}{2}} &= \epsilon \tilde{E}_l \cdot \tilde{E}_{l + 1} + \mu \hat{H}^2_{l + \frac{1}{2}} \\
|
|
\end{aligned}
|
|
$$
|
|
|
|
Rewriting the Poynting theorem in terms of the energy expressions,
|
|
|
|
$$
|
|
\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}
|
|
$$
|
|
|
|
This result is exact and should practically hold to within numerical precision. No time-
|
|
or spatial-averaging is necessary.
|
|
|
|
Note that each value of $J$ contributes to the energy twice (i.e. once per field update)
|
|
despite only causing the value of $E$ to change once (same for $M$ and $H$).
|
|
|
|
|
|
Sources
|
|
=============
|
|
|
|
It is often useful to excite the simulation with an arbitrary broadband pulse and then
|
|
extract the frequency-domain response by performing an on-the-fly Fourier transform
|
|
of the time-domain fields.
|
|
|
|
`accumulate_phasor` in `meanas.fdtd.phasor` performs the phase accumulation for one
|
|
or more target frequencies, while leaving source normalization and simulation-loop
|
|
policy to the caller. `temporal_phasor(...)` and `temporal_phasor_scale(...)`
|
|
apply the same Fourier sum to a scalar waveform, which is useful when a pulsed
|
|
source envelope must be normalized before being applied to a point source or
|
|
mode source. `real_injection_scale(...)` is the corresponding helper for the
|
|
common real-valued injection pattern `numpy.real(scale * waveform)`. Convenience
|
|
wrappers `accumulate_phasor_e`, `accumulate_phasor_h`, and `accumulate_phasor_j`
|
|
apply the usual Yee time offsets. `reconstruct_real(...)` and the corresponding
|
|
`reconstruct_real_e/h/j(...)` wrappers perform the inverse operation, converting
|
|
phasors back into real-valued field snapshots at explicit Yee-aligned times.
|
|
For scalar `omega`, the reconstruction helpers accept either a plain field
|
|
phasor or the batched `(1, *sample_shape)` form used by the accumulator helpers.
|
|
The helpers accumulate
|
|
|
|
$$ \Delta_t \sum_l w_l e^{-i \omega t_l} f_l $$
|
|
|
|
with caller-provided sample times and weights. In this codebase the matching
|
|
electric-current convention is typically `E -= dt * J / epsilon` in FDTD and
|
|
`-i \omega J` on the right-hand side of the FDFD wave equation.
|
|
|
|
For FDTD/FDFD equivalence, this phasor path is the primary comparison workflow.
|
|
It isolates the guided `+\omega` response that the frequency-domain solver
|
|
targets directly, regardless of whether the underlying FDTD run used real- or
|
|
complex-valued fields.
|
|
|
|
For exact pulsed FDTD/FDFD equivalence it is often simplest to keep the
|
|
injected source, fields, and CPML auxiliary state complex-valued. The
|
|
`real_injection_scale(...)` helper is instead for the more ordinary one-run
|
|
real-valued source path, where the intended positive-frequency waveform is
|
|
injected through `numpy.real(scale * waveform)` and any remaining negative-
|
|
frequency leakage is controlled by the pulse bandwidth and run length.
|
|
|
|
`reconstruct_real(...)` is for a different question: given a phasor, what late
|
|
real-valued field snapshot should it produce? That raw-snapshot comparison is
|
|
stricter and noisier because a monitor plane generally contains both the guided
|
|
field and the remaining orthogonal content,
|
|
|
|
$$ E_{\text{monitor}} = E_{\text{guided}} + E_{\text{residual}} . $$
|
|
|
|
Phasor/modal comparisons mostly validate the guided `+\omega` term. Raw
|
|
real-field comparisons expose both terms at once, so they should be treated as
|
|
secondary diagnostics rather than the main solver-equivalence benchmark.
|
|
|
|
The Ricker wavelet (normalized second derivative of a Gaussian) is commonly used for the pulse
|
|
shape. It can be written
|
|
|
|
$$ f_r(t) = (1 - \frac{1}{2} (\omega (t - \tau))^2) e^{-(\frac{\omega (t - \tau)}{2})^2} $$
|
|
|
|
with $\tau > \frac{2 * \pi}{\omega}$ as a minimum delay to avoid a discontinuity at
|
|
t=0 (assuming the source is off for t<0 this gives $\sim 10^{-3}$ error at t=0).
|
|
|
|
|
|
|
|
Boundary conditions
|
|
===================
|
|
|
|
`meanas.fdtd` exposes two boundary-related building blocks:
|
|
|
|
- `conducting_boundary(...)` for simple perfect-electric-conductor style field
|
|
clamping at one face of the domain.
|
|
- `cpml_params(...)` and `updates_with_cpml(...)` for convolutional perfectly
|
|
matched layers (CPMLs) attached to one or more faces of the Yee grid.
|
|
|
|
`updates_with_cpml(...)` accepts a three-by-two table of CPML parameter blocks:
|
|
|
|
```
|
|
cpml_params[axis][polarity_index]
|
|
```
|
|
|
|
where `axis` is `0`, `1`, or `2` and `polarity_index` corresponds to `(-1, +1)`.
|
|
Passing `None` for one entry disables CPML on that face while leaving the other
|
|
directions unchanged. This is how mixed boundary setups such as "absorbing in x,
|
|
periodic in y/z" are expressed.
|
|
|
|
When comparing an FDTD run against an FDFD solve, use the same stretched
|
|
coordinate system in both places:
|
|
|
|
1. Build the FDTD update with the desired CPML parameters.
|
|
2. Stretch the FDFD `dxes` with the matching SCPML transform.
|
|
3. Compare the extracted phasor against the FDFD field or residual on those
|
|
stretched `dxes`.
|
|
|
|
The electric-current sign convention used throughout the examples and tests is
|
|
|
|
$$
|
|
E \leftarrow E - \Delta_t J / \epsilon
|
|
$$
|
|
|
|
which matches the FDFD right-hand side
|
|
|
|
$$
|
|
-i \omega J.
|
|
$$
|
|
"""
|
|
|
|
from .base import (
|
|
maxwell_e as maxwell_e,
|
|
maxwell_h as maxwell_h,
|
|
)
|
|
from .pml import (
|
|
cpml_params as cpml_params,
|
|
updates_with_cpml as updates_with_cpml,
|
|
)
|
|
from .energy import (
|
|
poynting as poynting,
|
|
poynting_divergence as poynting_divergence,
|
|
energy_hstep as energy_hstep,
|
|
energy_estep as energy_estep,
|
|
delta_energy_h2e as delta_energy_h2e,
|
|
delta_energy_j as delta_energy_j,
|
|
)
|
|
from .boundaries import (
|
|
conducting_boundary as conducting_boundary,
|
|
)
|
|
from .phasor import (
|
|
accumulate_phasor as accumulate_phasor,
|
|
accumulate_phasor_e as accumulate_phasor_e,
|
|
accumulate_phasor_h as accumulate_phasor_h,
|
|
accumulate_phasor_j as accumulate_phasor_j,
|
|
real_injection_scale as real_injection_scale,
|
|
reconstruct_real as reconstruct_real,
|
|
reconstruct_real_e as reconstruct_real_e,
|
|
reconstruct_real_h as reconstruct_real_h,
|
|
reconstruct_real_j as reconstruct_real_j,
|
|
temporal_phasor as temporal_phasor,
|
|
temporal_phasor_scale as temporal_phasor_scale,
|
|
)
|