From f9d90378b4c336b3de9e88e6c6cce2dd9103a9af Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Wed, 8 Jan 2020 00:51:56 -0800 Subject: [PATCH] more documentation --- meanas/fdmath/__init__.py | 137 +++++++++++++++++++++++++++++++++++++- meanas/fdmath/types.py | 2 +- meanas/fdtd/__init__.py | 30 +++++++++ meanas/fdtd/base.py | 89 ++++++++++++++++++++++++- meanas/fdtd/boundaries.py | 2 + meanas/fdtd/energy.py | 7 ++ meanas/fdtd/pml.py | 3 + 7 files changed, 266 insertions(+), 4 deletions(-) diff --git a/meanas/fdmath/__init__.py b/meanas/fdmath/__init__.py index 0803099..76bd908 100644 --- a/meanas/fdmath/__init__.py +++ b/meanas/fdmath/__init__.py @@ -392,13 +392,108 @@ $$ -\\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}} \\\\ - \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l, \\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}} \\end{align*} $$ +Frequency domain +---------------- + +We can substitute in a time-harmonic fields + +$$ + \\begin{align*} + \\tilde{E}_\\vec{r} &= \\tilde{E}_\\vec{r} e^{-\\imath \\omega l \\Delta_t} \\\\ + \\tilde{J}_\\vec{r} &= \\tilde{J}_\\vec{r} e^{-\\imath \\omega (l - \\frac{1}{2}) \\Delta_t} + \\end{align*} +$$ + +resulting in + +$$ + \\begin{align*} + \\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{align*} +$$ + +This gives the frequency-domain wave equation, + +$$ + \\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} +$$ + + +Plane waves and Dispersion relation +------------------------------------ + +With uniform material distribution and no sources + +$$ + \\begin{align*} + \\mu_{\\vec{r} + \\frac{1}{2}} &= \\mu \\\\ + \\epsilon_\\vec{r} &= \\epsilon \\\\ + \\tilde{J}_\\vec{r} &= 0 \\\\ + \\end{align*} +$$ + +the frequency domain wave equation simplifies to + +$$ \\hat{\\nabla} \\times \\tilde{\\nabla} \\times \\tilde{E}_\\vec{r} - \\Omega^2 \\epsilon \\mu \\tilde{E}_\\vec{r} = 0 $$ + +Since \\( \\hat{\\nabla} \\cdot \\tilde{E}_\\vec{r} = 0 \\), we can simplify + +$$ + \\begin{align*} + \\hat{\\nabla} \\times \\tilde{\\nabla} \\times \\tilde{E}_\\vec{r} + &= \\tilde{\\nabla}(\\hat{\\nabla} \\cdot \\tilde{E}_\\vec{r}) - \\hat{\\nabla} \\cdot \\tilde{\\nabla} \\tilde{E}_\\vec{r} \\\\ + &= - \\hat{\\nabla} \\cdot \\tilde{\\nabla} \\tilde{E}_\\vec{r} \\\\ + &= - \\tilde{\\nabla}^2 \\tilde{E}_\\vec{r} + \\end{align*} +$$ + +and we get + +$$ \\tilde{\\nabla}^2 \\tilde{E}_\\vec{r} + \\Omega^2 \\epsilon \\mu \\tilde{E}_\\vec{r} = 0 $$ + +We can convert this to three scalar-wave equations of the form + +$$ (\\tilde{\\nabla}^2 + K^2) \\phi_\\vec{r} = 0 $$ + +with \\( K^2 = \\Omega^2 \\mu \\epsilon \\). Now we let + +$$ \\phi_\\vec{r} = A e^{\\imath (k_x m \\Delta_x + k_y n \\Delta_y + k_z p \\Delta_z)} $$ + +resulting in + +$$ + \\begin{align*} + \\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}\\\\ + K_x &= 2 \\sin(k_x \\Delta_x / 2) / \\Delta_x \\\\ + \\end{align*} +$$ + +with similar expressions for the y and z dimnsions (and \\( K_y, K_z \\)). + +This implies + +$$ + \\tilde{\\nabla}^2 = -(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 +$$ + +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}} \\). + Grid description ================ @@ -598,6 +693,46 @@ H-field back-vectors: where `dx_e[0]` is the x-width of the `m=0` cells, as used when calculating dE/dx, and `dy_h[0]` is the y-width of the `n=0` cells, as used when calculating dH/dy, etc. + +Permittivity and Permeability +============================= + +Since each vector component of E and H is defined in a different location and represents +a different volume, the value of the spatially-discrete `epsilon` and `mu` can also be +different for all three field components, even when representing a simple planar interface +between two isotropic materials. + +As a result, `epsilon` and `mu` are taken to have the same dimensions as the field, and +composed of the three diagonal tensor components: + + [equations: epsilon_and_mu] + epsilon = [epsilon_xx, epsilon_yy, epsilon_zz] + mu = [mu_xx, mu_yy, mu_zz] + +or + +$$ + \\epsilon = \\begin{bmatrix} \\epsilon_{xx} & 0 & 0 \\\\ + 0 & \\epsilon_{yy} & 0 \\\\ + 0 & 0 & \\epsilon_{zz} \\end{bmatrix} +$$ +$$ + \\mu = \\begin{bmatrix} \\mu_{xx} & 0 & 0 \\\\ + 0 & \\mu_{yy} & 0 \\\\ + 0 & 0 & \\mu_{zz} \\end{bmatrix} +$$ + +where the off-diagonal terms (e.g. `epsilon_xy`) are assumed to be zero. + +High-accuracy volumetric integration of shapes on multiple grids can be performed +by the [gridlock](https://mpxd.net/code/jan/gridlock) module. + +The values of the vacuum permittivity and permability effectively become scaling +factors that appear in several locations (e.g. between the E and H fields). In +order to limit floating-point inaccuracy and simplify calculations, they are often +set to 1 and relative permittivities and permeabilities are used in their places; +the true values can be multiplied back in after the simulation is complete if non- +normalized results are needed. """ from .types import fdfield_t, vfdfield_t, dx_lists_t, fdfield_updater_t diff --git a/meanas/fdmath/types.py b/meanas/fdmath/types.py index f8719d1..ea78f8f 100644 --- a/meanas/fdmath/types.py +++ b/meanas/fdmath/types.py @@ -37,7 +37,7 @@ dx_lists_t = List[List[numpy.ndarray]] ''' -fdfield_updater_t = Callable[[fdfield_t], fdfield_t] +fdfield_updater_t = Callable[..., fdfield_t] ''' Convenience type for functions which take and return an fdfield_t ''' diff --git a/meanas/fdtd/__init__.py b/meanas/fdtd/__init__.py index 02cfc05..2d1588e 100644 --- a/meanas/fdtd/__init__.py +++ b/meanas/fdtd/__init__.py @@ -1,5 +1,35 @@ """ Utilities for running finite-difference time-domain (FDTD) simulations + + +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 +=============== +# TODO + +Energy conservation +=================== +# TODO + +Boundary conditions +=================== +# TODO notes about boundaries / PMLs """ from .base import maxwell_e, maxwell_h diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py index fa478ba..fd21642 100644 --- a/meanas/fdtd/base.py +++ b/meanas/fdtd/base.py @@ -1,5 +1,7 @@ """ Basic FDTD field updates + + """ from typing import List, Callable, Tuple, Dict import numpy @@ -12,12 +14,51 @@ __author__ = 'Jan Petykiewicz' def maxwell_e(dt: float, dxes: dx_lists_t = None) -> fdfield_updater_t: + """ + Build a function which performs a portion the time-domain E-field update, + + E += curl_back(H[t]) / epsilon + + The full update should be + + E += (curl_back(H[t]) + J) / epsilon + + which requires an additional step of `E += J / epsilon` which is not performed + by the generated function. + + See `meanas.fdmath` for descriptions of + + - This update step: "Maxwell's equations" section + - `dxes`: "Datastructure: dx_lists_t" section + - `epsilon`: "Permittivity and Permeability" section + + Also see the "Timestep" section of `meanas.fdtd` for a discussion of + the `dt` parameter. + + Args: + dt: Timestep. See `meanas.fdtd` for details. + dxes: Grid description; see `meanas.fdmath`. + + Returns: + Function `f(E_old, H_old, epsilon) -> E_new`. + """ if dxes is not None: curl_h_fun = curl_back(dxes[1]) else: curl_h_fun = curl_back() def me_fun(e: fdfield_t, h: fdfield_t, epsilon: fdfield_t): + """ + Update the E-field. + + Args: + e: E-field at time t=0 + h: H-field at time t=0.5 + epsilon: Dielectric constant distribution. + + Returns: + E-field at time t=1 + """ e += dt * curl_h_fun(h) / epsilon return e @@ -25,13 +66,57 @@ def maxwell_e(dt: float, dxes: dx_lists_t = None) -> fdfield_updater_t: def maxwell_h(dt: float, dxes: dx_lists_t = None) -> fdfield_updater_t: + """ + Build a function which performs part of the time-domain H-field update, + + H -= curl_forward(E[t]) / mu + + The full update should be + + H -= (curl_forward(E[t]) - M) / mu + + 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`. + + See `meanas.fdmath` for descriptions of + + - This update step: "Maxwell's equations" section + - `dxes`: "Datastructure: dx_lists_t" section + - `mu`: "Permittivity and Permeability" section + + Also see the "Timestep" section of `meanas.fdtd` for a discussion of + the `dt` parameter. + + Args: + dt: Timestep. See `meanas.fdtd` for details. + dxes: Grid description; see `meanas.fdmath`. + + Returns: + Function `f(E_old, H_old, epsilon) -> E_new`. + """ if dxes is not None: curl_e_fun = curl_forward(dxes[0]) else: curl_e_fun = curl_forward() - def mh_fun(e: fdfield_t, h: fdfield_t): - h -= dt * curl_e_fun(e) + def mh_fun(e: fdfield_t, h: fdfield_t, mu: fdfield_t = None): + """ + Update the H-field. + + Args: + e: E-field at time t=1 + h: H-field at time t=0.5 + mu: Magnetic permeability. Default is 1 everywhere. + + Returns: + H-field at time t=1.5 + """ + if mu is not None: + h -= dt * curl_e_fun(e) / mu + else: + h -= dt * curl_e_fun(e) + return h return mh_fun diff --git a/meanas/fdtd/boundaries.py b/meanas/fdtd/boundaries.py index 0f0b1a6..10c966e 100644 --- a/meanas/fdtd/boundaries.py +++ b/meanas/fdtd/boundaries.py @@ -1,5 +1,7 @@ """ Boundary conditions + +#TODO conducting boundary documentation """ from typing import List, Callable, Tuple, Dict diff --git a/meanas/fdtd/energy.py b/meanas/fdtd/energy.py index 41268b7..dc5848a 100644 --- a/meanas/fdtd/energy.py +++ b/meanas/fdtd/energy.py @@ -5,10 +5,14 @@ import numpy from ..fdmath import dx_lists_t, fdfield_t, fdfield_updater_t from ..fdmath.functional import deriv_back, deriv_forward + def poynting(e: fdfield_t, h: fdfield_t, dxes: dx_lists_t = None, ) -> fdfield_t: + """ + Calculate the poynting vector + """ if dxes is None: dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2)) @@ -32,6 +36,9 @@ def poynting_divergence(s: fdfield_t = None, h: fdfield_t = None, dxes: dx_lists_t = None, ) -> fdfield_t: + """ + Calculate the divergence of the poynting vector + """ if s is None: s = poynting(e, h, dxes=dxes) diff --git a/meanas/fdtd/pml.py b/meanas/fdtd/pml.py index af15a5a..2c4ae1e 100644 --- a/meanas/fdtd/pml.py +++ b/meanas/fdtd/pml.py @@ -1,6 +1,9 @@ """ PML implementations +#TODO discussion of PMLs +#TODO cpml documentation + """ # TODO retest pmls!