From 34075a74560b19f665eaf60a6d7a4a830a03ae15 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 28 Oct 2026 17:49:32 -0700 Subject: [PATCH] [docs] correct Yee calculus, source polarity, and phasor conventions --- meanas/fdfd/__init__.py | 42 ++++++++++++++--------------- meanas/fdmath/__init__.py | 57 +++++++++++++++++++++------------------ meanas/fdtd/base.py | 12 ++++----- 3 files changed, 58 insertions(+), 53 deletions(-) diff --git a/meanas/fdfd/__init__.py b/meanas/fdfd/__init__.py index e94adc3..240d711 100644 --- a/meanas/fdfd/__init__.py +++ b/meanas/fdfd/__init__.py @@ -19,38 +19,39 @@ Submodules: ================================================================ -From the "Frequency domain" section of `meanas.fdmath`, we have +Following the "Frequency domain" section of `meanas.fdmath`, we use +$e^{+\imath \omega t}$ time dependence. Phasors are referenced to the +actual sample times of the staggered fields and sources: $$ \begin{aligned} - \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\ - \tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \tilde{H}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\ - \tilde{J}_{l, \vec{r}} &= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\ - \tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \tilde{M}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega l \Delta_t} \\ - \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}}) - -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} &= -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\ + \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{+\imath \omega l \Delta_t} \\ + \hat{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \hat{H}_{\vec{r} + \frac{1}{2}} e^{+\imath \omega (l - \frac{1}{2}) \Delta_t} \\ + \tilde{J}_{l - \frac{1}{2}, \vec{r}} &= \tilde{J}_{\vec{r}} e^{+\imath \omega (l - \frac{1}{2}) \Delta_t} \\ + \hat{M}_{l, \vec{r} + \frac{1}{2}} &= \hat{M}_{\vec{r} + \frac{1}{2}} e^{+\imath \omega l \Delta_t} \\ \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t \end{aligned} $$ -resulting in +The difference multipliers for an integer-indexed sequence are $$ \begin{aligned} - \tilde{\partial}_t &\Rightarrow -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\ - \hat{\partial}_t &\Rightarrow -\imath \Omega e^{ \imath \omega \Delta_t / 2}\\ + \tilde{\partial}_t &\Rightarrow \imath \Omega e^{ \imath \omega \Delta_t / 2}\\ + \hat{\partial}_t &\Rightarrow \imath \Omega e^{-\imath \omega \Delta_t / 2}\\ \end{aligned} $$ -Maxwell's equations are then +The half-step factors cancel when the staggered sample times are included. +Maxwell's equations for the phasors are then $$ \begin{aligned} \tilde{\nabla} \times \tilde{E}_{\vec{r}} &= - \imath \Omega e^{-\imath \omega \Delta_t / 2} \hat{B}_{\vec{r} + \frac{1}{2}} + -\imath \Omega \hat{B}_{\vec{r} + \frac{1}{2}} - \hat{M}_{\vec{r} + \frac{1}{2}} \\ \hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &= - -\imath \Omega e^{ \imath \omega \Delta_t / 2} \tilde{D}_{\vec{r}} + \imath \Omega \tilde{D}_{\vec{r}} + \tilde{J}_{\vec{r}} \\ \tilde{\nabla} \cdot \hat{B}_{\vec{r} + \frac{1}{2}} &= 0 \\ \hat{\nabla} \cdot \tilde{D}_{\vec{r}} &= \rho_{\vec{r}} @@ -61,13 +62,9 @@ With $\Delta_t \to 0$, this simplifies to $$ \begin{aligned} - \tilde{E}_{l, \vec{r}} &\to \tilde{E}_{\vec{r}} \\ - \tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &\to \tilde{H}_{\vec{r} + \frac{1}{2}} \\ - \tilde{J}_{l, \vec{r}} &\to \tilde{J}_{\vec{r}} \\ - \tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &\to \tilde{M}_{\vec{r} + \frac{1}{2}} \\ \Omega &\to \omega \\ - \tilde{\partial}_t &\to -\imath \omega \\ - \hat{\partial}_t &\to -\imath \omega \\ + \tilde{\partial}_t &\to \imath \omega \\ + \hat{\partial}_t &\to \imath \omega \\ \end{aligned} $$ @@ -76,14 +73,17 @@ and then $$ \begin{aligned} \tilde{\nabla} \times \tilde{E}_{\vec{r}} &= - \imath \omega \hat{B}_{\vec{r} + \frac{1}{2}} + -\imath \omega \hat{B}_{\vec{r} + \frac{1}{2}} - \hat{M}_{\vec{r} + \frac{1}{2}} \\ \hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &= - -\imath \omega \tilde{D}_{\vec{r}} + \imath \omega \tilde{D}_{\vec{r}} + \tilde{J}_{\vec{r}} \\ \end{aligned} $$ +For zero magnetic current, eliminating $H$ gives the electric wave equation +used by `operators.e_full` and `functional.e_full`: + $$ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}}) -\omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \omega \tilde{J}_{\vec{r}} \\ diff --git a/meanas/fdmath/__init__.py b/meanas/fdmath/__init__.py index 857cf18..ce9827b 100644 --- a/meanas/fdmath/__init__.py +++ b/meanas/fdmath/__init__.py @@ -116,9 +116,9 @@ Expanding to three dimensions, we can define two gradients $$ [\tilde{\nabla} f]_{m,n,p} = \vec{x} [\tilde{\partial}_x f]_{m + \frac{1}{2},n,p} + \vec{y} [\tilde{\partial}_y f]_{m,n + \frac{1}{2},p} + \vec{z} [\tilde{\partial}_z f]_{m,n,p + \frac{1}{2}} $$ - $$ [\hat{\nabla} f]_{m,n,p} = \vec{x} [\hat{\partial}_x f]_{m + \frac{1}{2},n,p} + - \vec{y} [\hat{\partial}_y f]_{m,n + \frac{1}{2},p} + - \vec{z} [\hat{\partial}_z f]_{m,n,p + \frac{1}{2}} $$ + $$ [\hat{\nabla} f]_{m,n,p} = \vec{x} [\hat{\partial}_x f]_{m - \frac{1}{2},n,p} + + \vec{y} [\hat{\partial}_y f]_{m,n - \frac{1}{2},p} + + \vec{z} [\hat{\partial}_z f]_{m,n,p - \frac{1}{2}} $$ or @@ -232,7 +232,7 @@ The two curls are then [\tilde{\nabla} \times \tilde{g}]_{m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2}} &= \vec{x} (\tilde{\partial}_y g^z_{m,n,p + \frac{1}{2}} - \tilde{\partial}_z g^y_{m,n + \frac{1}{2},p}) \\ &+ \vec{y} (\tilde{\partial}_z g^x_{m + \frac{1}{2},n,p} - \tilde{\partial}_x g^z_{m,n,p + \frac{1}{2}}) \\ - &+ \vec{z} (\tilde{\partial}_x g^y_{m,n + \frac{1}{2},p} - \tilde{\partial}_y g^z_{m + \frac{1}{2},n,p}) + &+ \vec{z} (\tilde{\partial}_x g^y_{m,n + \frac{1}{2},p} - \tilde{\partial}_y g^x_{m + \frac{1}{2},n,p}) \end{aligned} $$ and @@ -314,7 +314,7 @@ The time derivatives can be expanded to form the update equations: [code: Maxwell's equations updates] H[i, j, k] -= dt * (curl_forward(E)[i, j, k] + M[t, i, j, k]) / mu[i, j, k] - E[i, j, k] += dt * (curl_back( H)[i, j, k] + J[t, i, j, k]) / epsilon[i, j, k] + E[i, j, k] += dt * (curl_back( H)[i, j, k] - J[t, i, j, k]) / epsilon[i, j, k] Note that the E-field fore-vector and H-field back-vector are offset by a half-cell, resulting in distinct locations for all six E- and H-field components: @@ -384,7 +384,7 @@ $$ \begin{aligned} \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &= -\tilde{\partial}_t \hat{B}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}} - - \hat{M}_{l-1, \vec{r} + \frac{1}{2}} \\ + - \hat{M}_{l, \vec{r} + \frac{1}{2}} \\ \mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &= -\tilde{\partial}_t \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}} \\ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &= @@ -392,10 +392,10 @@ $$ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &= -\tilde{\partial}_t \hat{\nabla} \times \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}} \\ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &= - -\tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \tilde{E}_{l, \vec{r}} + \hat{\partial}_t \tilde{J}_{l-\frac{1}{2},\vec{r}} \\ + -\tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \tilde{E}_{l, \vec{r}} - \tilde{\partial}_t \tilde{J}_{l-\frac{1}{2},\vec{r}} \\ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) + \tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \cdot \tilde{E}_{l, \vec{r}} - &= \tilde{\partial}_t \tilde{J}_{l - \frac{1}{2}, \vec{r}} + &= -\tilde{\partial}_t \tilde{J}_{l - \frac{1}{2}, \vec{r}} \end{aligned} $$ @@ -403,37 +403,42 @@ $$ Frequency domain ---------------- -We can substitute in a time-harmonic fields +We use the $e^{+\imath \omega t}$ convention, matching the FDFD operators. +Substitute time-harmonic fields with phases referenced to their actual sample times: $$ \begin{aligned} - \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\ - \tilde{J}_{l, \vec{r}} &= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} + \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{+\imath \omega l \Delta_t} \\ + \tilde{J}_{l-\frac{1}{2}, \vec{r}} &= \tilde{J}_{\vec{r}} e^{+\imath \omega (l - \frac{1}{2}) \Delta_t} \end{aligned} $$ -resulting in +For an integer-indexed sequence proportional to $e^{+\imath \omega l \Delta_t}$, +the forward and backward difference multipliers are $$ \begin{aligned} - \tilde{\partial}_t &\Rightarrow (e^{ \imath \omega \Delta_t} - 1) / \Delta_t = \frac{-2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{-\imath \omega \Delta_t / 2} = -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\ - \hat{\partial}_t &\Rightarrow (1 - e^{-\imath \omega \Delta_t}) / \Delta_t = \frac{-2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{ \imath \omega \Delta_t / 2} = -\imath \Omega e^{ \imath \omega \Delta_t / 2}\\ + \tilde{\partial}_t &\Rightarrow (e^{ \imath \omega \Delta_t} - 1) / \Delta_t = \frac{2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{ \imath \omega \Delta_t / 2} = \imath \Omega e^{ \imath \omega \Delta_t / 2}\\ + \hat{\partial}_t &\Rightarrow (1 - e^{-\imath \omega \Delta_t}) / \Delta_t = \frac{2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{-\imath \omega \Delta_t / 2} = \imath \Omega e^{-\imath \omega \Delta_t / 2}\\ \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t \end{aligned} $$ -This gives the frequency-domain wave equation, +The half-step factors account for the shift in the derivative's sample time. +In particular, $\tilde{\partial}_t J_{l-\frac{1}{2}} = \imath \Omega J e^{\imath \omega l \Delta_t}$. +Thus, when the staggered fields are expressed using the phasors above, the +frequency-domain wave equation has no extra half-step phase: $$ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}}) - -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\ + -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \Omega \tilde{J}_{\vec{r}} \\ $$ Plane waves and Dispersion relation ------------------------------------ -With uniform material distribution and no sources +With uniform grid spacings, a uniform material distribution, and no sources $$ \begin{aligned} @@ -474,31 +479,31 @@ resulting in $$ \begin{aligned} - \tilde{\partial}_x &\Rightarrow (e^{ \imath k_x \Delta_x} - 1) / \Delta_t = \frac{-2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{ \imath k_x \Delta_x / 2} = \imath K_x e^{ \imath k_x \Delta_x / 2}\\ - \hat{\partial}_x &\Rightarrow (1 - e^{-\imath k_x \Delta_x}) / \Delta_t = \frac{-2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{-\imath k_x \Delta_x / 2} = \imath K_x e^{-\imath k_x \Delta_x / 2}\\ + \tilde{\partial}_x &\Rightarrow (e^{ \imath k_x \Delta_x} - 1) / \Delta_x = \frac{2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{ \imath k_x \Delta_x / 2} = \imath K_x e^{ \imath k_x \Delta_x / 2}\\ + \hat{\partial}_x &\Rightarrow (1 - e^{-\imath k_x \Delta_x}) / \Delta_x = \frac{2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{-\imath k_x \Delta_x / 2} = \imath K_x e^{-\imath k_x \Delta_x / 2}\\ K_x &= 2 \sin(k_x \Delta_x / 2) / \Delta_x \\ \end{aligned} $$ -with similar expressions for the y and z dimnsions (and $K_y, K_z$). +with similar expressions for the y and z dimensions (and $K_y, K_z$). This implies $$ - \tilde{\nabla}^2 = -(K_x^2 + K_y^2 + K_z^2) \phi_{\vec{r}} \\ + \tilde{\nabla}^2 \phi_{\vec{r}} = -(K_x^2 + K_y^2 + K_z^2) \phi_{\vec{r}} \\ K_x^2 + K_y^2 + K_z^2 = \Omega^2 \mu \epsilon = \Omega^2 / c^2 $$ -where $c = \sqrt{\mu \epsilon}$. +where $c = 1 / \sqrt{\mu \epsilon}$. -Assuming real $(k_x, k_y, k_z), \omega$ will be real only if +To keep $\omega$ real for every real grid wavevector $(k_x, k_y, k_z)$, +we use the conservative stability bound $$ 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}) $$ If $\Delta_x = \Delta_y = \Delta_z$, this simplifies to $c \Delta_t < \Delta_x / \sqrt{3}$. -This last form can be interpreted as enforcing causality; the distance that light -travels in one timestep (i.e., $c \Delta_t$) must be less than the diagonal -of the smallest cell ( $\Delta_x / \sqrt{3}$ when on a uniform cubic grid). +This limits the distance that light travels in one timestep to less than +$\Delta_x / \sqrt{3}$ on a uniform cubic grid. Grid description diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py index 480ed87..769f362 100644 --- a/meanas/fdtd/base.py +++ b/meanas/fdtd/base.py @@ -17,13 +17,13 @@ def maxwell_e( """ Build a function which performs a portion the time-domain E-field update, - E += curl_back(H[t]) / epsilon + E += dt * curl_back(H[t]) / epsilon The full update should be - E += (curl_back(H[t]) + J) / epsilon + E += dt * (curl_back(H[t]) - J) / epsilon - which requires an additional step of `E += J / epsilon` which is not performed + which requires an additional step of `E -= dt * J / epsilon` which is not performed by the generated function. See `meanas.fdmath` for descriptions of @@ -72,13 +72,13 @@ def maxwell_h( """ Build a function which performs part of the time-domain H-field update, - H -= curl_forward(E[t]) / mu + H -= dt * curl_forward(E[t]) / mu The full update should be - H -= (curl_forward(E[t]) + M) / mu + H -= dt * (curl_forward(E[t]) + M) / mu - which requires an additional step of `H -= M / mu` which is not performed + which requires an additional step of `H -= dt * M / mu` which is not performed by the generated function; this step can be omitted if there is no magnetic current `M`.