From 034f79eae6c335679fb1d96657d95f44ca8d8e02 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 12 Jan 2020 22:50:01 -0800 Subject: [PATCH] Poynting vector doc updates --- meanas/fdmath/__init__.py | 21 ++++--- meanas/fdtd/__init__.py | 122 ++++++++++++++++++++++++++++++++++++-- meanas/fdtd/base.py | 4 +- 3 files changed, 131 insertions(+), 16 deletions(-) diff --git a/meanas/fdmath/__init__.py b/meanas/fdmath/__init__.py index 76bd908..cb47c79 100644 --- a/meanas/fdmath/__init__.py +++ b/meanas/fdmath/__init__.py @@ -288,9 +288,9 @@ If we discretize both space (m,n,p) and time (l), Maxwell's equations become $$ \\begin{align*} \\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{\\nabla} \\times \\hat{H}_{l,\\vec{r} + \\frac{1}{2}} &= \\hat{\\partial}_t \\tilde{D}_{l, \\vec{r}} - + \\tilde{J}_{l-\\frac{1}{2},\\vec{r}} \\\\ + - \\hat{M}_{l, \\vec{r} + \\frac{1}{2}} \\\\ + \\hat{\\nabla} \\times \\hat{H}_{l-\\frac{1}{2},\\vec{r} + \\frac{1}{2}} &= \\hat{\\partial}_t \\tilde{D}_{l, \\vec{r}} + + \\tilde{J}_{l-\\frac{1}{2},\\vec{r}} \\\\ \\tilde{\\nabla} \\cdot \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &= 0 \\\\ \\hat{\\nabla} \\cdot \\tilde{D}_{l,\\vec{r}} &= \\rho_{l,\\vec{r}} \\end{align*} $$ @@ -311,9 +311,9 @@ and \\( \\epsilon \\) and \\( \\mu \\) are the dielectric permittivity and magne The above is Yee's algorithm, written in a form analogous to Maxwell's equations. The time derivatives can be expanded to form the update equations: - [code: Maxwell's equations] - H[i, j, k] -= (curl_forward(E[t])[i, j, k] - M[t, i, j, k]) / mu[i, j, k] - E[i, j, k] += (curl_back( H[t])[i, j, k] + J[t, i, j, k]) / epsilon[i, j, k] + [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] 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: @@ -383,7 +383,7 @@ $$ \\begin{align*} \\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-1, \\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}}) &= @@ -488,11 +488,16 @@ $$ K_x^2 + K_y^2 + K_z^2 = \\Omega^2 \\mu \\epsilon = \\Omega^2 / c^2 $$ +where \\( c = \\sqrt{\\mu \\epsilon} \\). + Assuming real \\( (k_x, k_y, k_z), \\omega \\) will be real only if $$ 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 < \\frac{\\Delta_x}{\\sqrt{3}} \\). +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). Grid description diff --git a/meanas/fdtd/__init__.py b/meanas/fdtd/__init__.py index 2d1588e..3d498a9 100644 --- a/meanas/fdtd/__init__.py +++ b/meanas/fdtd/__init__.py @@ -1,6 +1,9 @@ """ Utilities for running finite-difference time-domain (FDTD) simulations +See the discussion of `Maxwell's Equations` in `meanas.fdmath` for basic +mathematical background. + Timestep ======== @@ -19,13 +22,120 @@ Based on this, we can set The `dx_min`, `dy_min`, `dz_min` should be the minimum value across both the base and derived grids. -Poynting Vector -=============== -# TODO +Poynting Vector and Energy Conservation +======================================= + +Let + +$$ \\begin{align*} + \\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{align*} +$$ + +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{align*} + \\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-1, \\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-1} - \\tilde{E}_l \\cdot \\tilde{J}_{l'} \\\\ + \\end{align*} +$$ + +where in the last line the spatial subscripts have been dropped to emphasize +the time subscripts \\( l, l' \\), i.e. + +$$ + \\begin{align*} + \\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{align*} +$$ + +etc. +For \\( l' = l + \\frac{1}{2} \\) we get + +$$ + \\begin{align*} + \\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{align*} +$$ + +and for \\( l' = l - \\frac{1}{2} \\), + +$$ + \\begin{align*} + \\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{align*} +$$ + +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{align*} + 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{align*} +$$ + +Rewriting the Poynting theorem in terms of the energy expressions, + +$$ + \\begin{align*} + (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{align*} +$$ + +This result is exact an 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 \\)). -Energy conservation -=================== -# TODO Boundary conditions =================== diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py index fd21642..e16a53f 100644 --- a/meanas/fdtd/base.py +++ b/meanas/fdtd/base.py @@ -73,9 +73,9 @@ def maxwell_h(dt: float, dxes: dx_lists_t = None) -> fdfield_updater_t: The full update should be - H -= (curl_forward(E[t]) - M) / mu + H -= (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 -= M / mu` which is not performed by the generated function; this step can be omitted if there is no magnetic current `M`.