more documentation

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

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
 '''

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

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

meanas/fdtd/boundaries.py

@@ -1,5 +1,7 @@
 """
 Boundary conditions
+
+#TODO conducting boundary documentation
 """
 
 from typing import List, Callable, Tuple, Dict

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)
 

meanas/fdtd/pml.py

@@ -1,6 +1,9 @@
 """
 PML implementations
 
+#TODO discussion of PMLs
+#TODO cpml documentation
+
 """
 # TODO retest pmls!