jan/fdfd_tools
2
2
Fork 1
Code Issues Pull Requests Releases Wiki Activity
034f79eae6
fdfd_tools/meanas/fdmath/__init__.py

747 lines
36 KiB
Python
Raw Normal View History

documentation updates
2019-11-26 01:47:52 -08:00
"""
add headings and vector diagram
2019-11-30 01:24:16 -08:00
documentation updates
2019-11-26 01:47:52 -08:00
Basic discrete calculus for finite difference (fd) simulations.
add info about functions vs operators
2020-01-04 18:46:28 -08:00
Fields, Functions, and Operators
================================
Discrete fields are stored in one of two forms:
Update some type info
2020-01-06 00:08:08 -08:00
- The `fdfield_t` form is a multidimensional `numpy.ndarray`
add info about functions vs operators
2020-01-04 18:46:28 -08:00
+ For a scalar field, this is just `U[m, n, p]`, where `m`, `n`, and `p` are
discrete indices referring to positions on the x, y, and z axes respectively.
+ For a vector field, the first index specifies which vector component is accessed:
`E[:, m, n, p] = [Ex[m, n, p], Ey[m, n, p], Ez[m, n, p]]`.
Update some type info
2020-01-06 00:08:08 -08:00
- The `vfdfield_t` form is simply a vectorzied (i.e. 1D) version of the `field_t`,
add info about functions vs operators
2020-01-04 18:46:28 -08:00
as obtained by `meanas.fdmath.vectorization.vec` (effectively just `numpy.ravel`)
Operators which act on fields also come in two forms:
+ Python functions, created by the functions in `meanas.fdmath.functional`.
Update some type info
2020-01-06 00:08:08 -08:00
The generated functions act on fields in the `fdfield_t` form.
add info about functions vs operators
2020-01-04 18:46:28 -08:00
+ Linear operators, usually 2D sparse matrices using `scipy.sparse`, created
by `meanas.fdmath.operators`. These operators act on vectorized fields in the
Update some type info
2020-01-06 00:08:08 -08:00
`vfdfield_t` form.
add info about functions vs operators
2020-01-04 18:46:28 -08:00
The operations performed should be equivalent: `functional.op(*args)(E)` should be
equivalent to `unvec(operators.op(*args) @ vec(E), E.shape[1:])`.
Generally speaking the `field_t` form is easier to work with, but can be harder or less
efficient to compose (e.g. it is easy to generate a single matrix by multiplying a
series of other matrices).
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
add headings and vector diagram
2019-11-30 01:24:16 -08:00
Discrete calculus
=================
documentation updates
2019-11-26 01:47:52 -08:00
This documentation and approach is roughly based on W.C. Chew's excellent
"Electromagnetic Theory on a Lattice" (doi:10.1063/1.355770),
which covers a superset of this material with similar notation and more detail.
more grid work
2019-12-09 21:28:26 -08:00
Scalar derivatives and cell shifts
----------------------------------
add headings and vector diagram
2019-11-30 01:24:16 -08:00
documentation updates
2019-11-26 01:47:52 -08:00
Define the discrete forward derivative as
lots more docs
2019-12-01 02:32:31 -08:00
$$ [\\tilde{\\partial}_x f ]_{m + \\frac{1}{2}} = \\frac{1}{\\Delta_{x, m}} (f_{m + 1} - f_m) $$
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
where \\( f \\) is a function defined at discrete locations on the x-axis (labeled using \\( m \\)).
The value at \\( m \\) occupies a length \\( \\Delta_{x, m} \\) along the x-axis. Note that \\( m \\)
is an index along the x-axis, _not_ necessarily an x-coordinate, since each length
\\( \\Delta_{x, m}, \\Delta_{x, m+1}, ...\\) is independently chosen.
If we treat `f` as a 1D array of values, with the `i`-th value `f[i]` taking up a length `dx[i]`
along the x-axis, the forward derivative is
deriv_forward(f)[i] = (f[i + 1] - f[i]) / dx[i]
documentation updates
2019-11-26 01:47:52 -08:00
Likewise, discrete reverse derivative is
lots more docs
2019-12-01 02:32:31 -08:00
$$ [\\hat{\\partial}_x f ]_{m - \\frac{1}{2}} = \\frac{1}{\\Delta_{x, m}} (f_{m} - f_{m - 1}) $$
add headings and vector diagram
2019-11-30 01:24:16 -08:00
or
documentation updates
2019-11-26 01:47:52 -08:00
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
deriv_back(f)[i] = (f[i] - f[i - 1]) / dx[i]
The derivatives' values are shifted by a half-cell relative to the original function, and
will have different cell widths if all the `dx[i]` ( \\( \\Delta_{x, m} \\) ) are not
identical:
[figure: derivatives and cell sizes]
dx0 dx1 dx2 dx3 cell sizes for function
----- ----- ----------- -----
______________________________
| | | |
f0 | f1 | f2 | f3 | function
_____|_____|___________|_____|
| | | |
| Df0 | Df1 | Df2 | Df3 forward derivative (periodic boundary)
__|_____|________|________|___
dx'3] dx'0 dx'1 dx'2 [dx'3 cell sizes for forward derivative
-- ----- -------- -------- ---
dx'0] dx'1 dx'2 dx'3 [dx'0 cell sizes for reverse derivative
______________________________
| | | |
| df1 | df2 | df3 | df0 reverse derivative (periodic boundary)
__|_____|________|________|___
Periodic boundaries are used here and elsewhere unless otherwise noted.
In the above figure,
`f0 =` \\(f_0\\), `f1 =` \\(f_1\\)
`Df0 =` \\([\\tilde{\\partial}f]_{0 + \\frac{1}{2}}\\)
`Df1 =` \\([\\tilde{\\partial}f]_{1 + \\frac{1}{2}}\\)
`df0 =` \\([\\hat{\\partial}f]_{0 - \\frac{1}{2}}\\)
etc.
The fractional subscript \\( m + \\frac{1}{2} \\) is used to indicate values defined
at shifted locations relative to the original \\( m \\), with corresponding lengths
$$ \\Delta_{x, m + \\frac{1}{2}} = \\frac{1}{2} * (\\Delta_{x, m} + \\Delta_{x, m + 1}) $$
Just as \\( m \\) is not itself an x-coordinate, neither is \\( m + \\frac{1}{2} \\);
carefully note the positions of the various cells in the above figure vs their labels.
more grid work
2019-12-09 21:28:26 -08:00
If the positions labeled with \\( m \\) are considered the "base" or "original" grid,
the positions labeled with \\( m + \\frac{1}{2} \\) are said to lie on a "dual" or
"derived" grid.
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
For the remainder of the `Discrete calculus` section, all figures will show
constant-length cells in order to focus on the vector derivatives themselves.
more grid work
2019-12-09 21:28:26 -08:00
See the `Grid description` section below for additional information on this topic
and generalization to three dimensions.
documentation updates
2019-11-26 01:47:52 -08:00
add headings and vector diagram
2019-11-30 01:24:16 -08:00
Gradients and fore-vectors
--------------------------
documentation updates
2019-11-26 01:47:52 -08:00
Expanding to three dimensions, we can define two gradients
lots more docs
2019-12-01 02:32:31 -08:00
$$ [\\tilde{\\nabla} f]_{m,n,p} = \\vec{x} [\\tilde{\\partial}_x f]_{m + \\frac{1}{2},n,p} +
documentation updates
2019-11-26 01:47:52 -08:00
\\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}} $$
lots more docs
2019-12-01 02:32:31 -08:00
or
[code: gradients]
grad_forward(f)[i,j,k] = [Dx_forward(f)[i, j, k],
Dy_forward(f)[i, j, k],
Dz_forward(f)[i, j, k]]
= [(f[i + 1, j, k] - f[i, j, k]) / dx[i],
(f[i, j + 1, k] - f[i, j, k]) / dy[i],
(f[i, j, k + 1] - f[i, j, k]) / dz[i]]
grad_back(f)[i,j,k] = [Dx_back(f)[i, j, k],
Dy_back(f)[i, j, k],
Dz_back(f)[i, j, k]]
= [(f[i, j, k] - f[i - 1, j, k]) / dx[i],
(f[i, j, k] - f[i, j - 1, k]) / dy[i],
(f[i, j, k] - f[i, j, k - 1]) / dz[i]]
documentation updates
2019-11-26 01:47:52 -08:00
The three derivatives in the gradient cause shifts in different
directions, so the x/y/z components of the resulting "vector" are defined
at different points: the x-component is shifted in the x-direction,
y in y, and z in z.
We call the resulting object a "fore-vector" or "back-vector", depending
on the direction of the shift. We write it as
$$ \\tilde{g}_{m,n,p} = \\vec{x} g^x_{m + \\frac{1}{2},n,p} +
\\vec{y} g^y_{m,n + \\frac{1}{2},p} +
\\vec{z} g^z_{m,n,p + \\frac{1}{2}} $$
$$ \\hat{g}_{m,n,p} = \\vec{x} g^x_{m - \\frac{1}{2},n,p} +
\\vec{y} g^y_{m,n - \\frac{1}{2},p} +
\\vec{z} g^z_{m,n,p - \\frac{1}{2}} $$
lots more docs
2019-12-01 02:32:31 -08:00
[figure: gradient / fore-vector]
(m, n+1, p+1) ______________ (m+1, n+1, p+1)
/: /|
/ : / |
/ : / |
(m, n, p+1)/_____________/ | The forward derivatives are defined
| : | | at the Dx, Dy, Dz points,
| :.........|...| but the forward-gradient fore-vector
more grid work
2019-12-09 21:28:26 -08:00
z y Dz / | / is the set of all three
|/_x | Dy | / and is said to be "located" at (m,n,p)
|/ |/
(m, n, p)|_____Dx______| (m+1, n, p)
add headings and vector diagram
2019-11-30 01:24:16 -08:00
Divergences
-----------
documentation updates
2019-11-26 01:47:52 -08:00
There are also two divergences,
$$ d_{n,m,p} = [\\tilde{\\nabla} \\cdot \\hat{g}]_{n,m,p}
= [\\tilde{\\partial}_x g^x]_{m,n,p} +
[\\tilde{\\partial}_y g^y]_{m,n,p} +
[\\tilde{\\partial}_z g^z]_{m,n,p} $$
$$ d_{n,m,p} = [\\hat{\\nabla} \\cdot \\tilde{g}]_{n,m,p}
= [\\hat{\\partial}_x g^x]_{m,n,p} +
[\\hat{\\partial}_y g^y]_{m,n,p} +
[\\hat{\\partial}_z g^z]_{m,n,p} $$
lots more docs
2019-12-01 02:32:31 -08:00
or
[code: divergences]
div_forward(g)[i,j,k] = Dx_forward(gx)[i, j, k] +
Dy_forward(gy)[i, j, k] +
Dz_forward(gz)[i, j, k]
= (gx[i + 1, j, k] - gx[i, j, k]) / dx[i] +
(gy[i, j + 1, k] - gy[i, j, k]) / dy[i] +
(gz[i, j, k + 1] - gz[i, j, k]) / dz[i]
div_back(g)[i,j,k] = Dx_back(gx)[i, j, k] +
Dy_back(gy)[i, j, k] +
Dz_back(gz)[i, j, k]
= (gx[i, j, k] - gx[i - 1, j, k]) / dx[i] +
(gy[i, j, k] - gy[i, j - 1, k]) / dy[i] +
(gz[i, j, k] - gz[i, j, k - 1]) / dz[i]
where `g = [gx, gy, gz]` is a fore- or back-vector field.
documentation updates
2019-11-26 01:47:52 -08:00
Since we applied the forward divergence to the back-vector (and vice-versa), the resulting scalar value
is defined at the back-vector's (fore-vectors) location \\( (m,n,p) \\) and not at the locations of its components
\\( (m \\pm \\frac{1}{2},n,p) \\) etc.
lots more docs
2019-12-01 02:32:31 -08:00
[figure: divergence]
^^
(m-1/2, n+1/2, p+1/2) _____||_______ (m+1/2, n+1/2, p+1/2)
/: || ,, /|
/ : || // / | The divergence at (m, n, p) (the center
/ : // / | of this cube) of a fore-vector field
(m-1/2, n-1/2, p+1/2)/_____________/ | is the sum of the outward-pointing
| : | | fore-vector components, which are
more grid work
2019-12-09 21:28:26 -08:00
z y <==|== :.........|.====> located at the face centers.
|/_x | / | /
| / // | / Note that in a nonuniform grid, each
|/ // || |/ dimension is normalized by the cell width.
(m-1/2, n-1/2, p-1/2)|____//_______| (m+1/2, n-1/2, p-1/2)
lots more docs
2019-12-01 02:32:31 -08:00
'' ||
VV
documentation updates
2019-11-26 01:47:52 -08:00
add headings and vector diagram
2019-11-30 01:24:16 -08:00
Curls
-----
documentation updates
2019-11-26 01:47:52 -08:00
The two curls are then
add headings and vector diagram
2019-11-30 01:24:16 -08:00
lots more docs
2019-12-01 02:32:31 -08:00
$$ \\begin{align*}
documentation updates
2019-11-26 01:47:52 -08:00
\\hat{h}_{m + \\frac{1}{2}, n + \\frac{1}{2}, p + \\frac{1}{2}} &= \\\\
[\\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}}) \\\\
lots more docs
2019-12-01 02:32:31 -08:00
&+ \\vec{z} (\\tilde{\\partial}_x g^y_{m,n + \\frac{1}{2},p} - \\tilde{\\partial}_y g^z_{m + \\frac{1}{2},n,p})
\\end{align*} $$
add headings and vector diagram
2019-11-30 01:24:16 -08:00
and
documentation updates
2019-11-26 01:47:52 -08:00
$$ \\tilde{h}_{m - \\frac{1}{2}, n - \\frac{1}{2}, p - \\frac{1}{2}} =
[\\hat{\\nabla} \\times \\hat{g}]_{m - \\frac{1}{2}, n - \\frac{1}{2}, p - \\frac{1}{2}} $$
where \\( \\hat{g} \\) and \\( \\tilde{g} \\) are located at \\((m,n,p)\\)
with components at \\( (m \\pm \\frac{1}{2},n,p) \\) etc.,
more doc fixes
2019-11-28 01:27:10 -08:00
while \\( \\hat{h} \\) and \\( \\tilde{h} \\) are located at \\((m \\pm \\frac{1}{2}, n \\pm \\frac{1}{2}, p \\pm \\frac{1}{2})\\)
documentation updates
2019-11-26 01:47:52 -08:00
with components at \\((m, n \\pm \\frac{1}{2}, p \\pm \\frac{1}{2})\\) etc.
lots more docs
2019-12-01 02:32:31 -08:00
[code: curls]
curl_forward(g)[i,j,k] = [Dy_forward(gz)[i, j, k] - Dz_forward(gy)[i, j, k],
Dz_forward(gx)[i, j, k] - Dx_forward(gz)[i, j, k],
Dx_forward(gy)[i, j, k] - Dy_forward(gx)[i, j, k]]
curl_back(g)[i,j,k] = [Dy_back(gz)[i, j, k] - Dz_back(gy)[i, j, k],
Dz_back(gx)[i, j, k] - Dx_back(gz)[i, j, k],
Dx_back(gy)[i, j, k] - Dy_back(gx)[i, j, k]]
For example, consider the forward curl, at (m, n, p), of a back-vector field `g`, defined
on a grid containing (m + 1/2, n + 1/2, p + 1/2).
The curl will be a fore-vector, so its z-component will be defined at (m, n, p + 1/2).
Take the nearest x- and y-components of `g` in the xy plane where the curl's z-component
is located; these are
[curl components]
(m, n + 1/2, p + 1/2) : x-component of back-vector at (m + 1/2, n + 1/2, p + 1/2)
(m + 1, n + 1/2, p + 1/2) : x-component of back-vector at (m + 3/2, n + 1/2, p + 1/2)
(m + 1/2, n , p + 1/2) : y-component of back-vector at (m + 1/2, n + 1/2, p + 1/2)
(m + 1/2, n + 1 , p + 1/2) : y-component of back-vector at (m + 1/2, n + 3/2, p + 1/2)
These four xy-components can be used to form a loop around the curl's z-component; its magnitude and sign
is set by their loop-oriented sum (i.e. two have their signs flipped to complete the loop).
[figure: z-component of curl]
more grid work
2019-12-09 21:28:26 -08:00
: |
z y : ^^ |
|/_x :....||.<.....| (m, n+1, p+1/2)
/ || /
| v || | ^
|/ |/
(m, n, p+1/2) |_____>______| (m+1, n, p+1/2)
lots more docs
2019-12-01 02:32:31 -08:00
Maxwell's Equations
===================
If we discretize both space (m,n,p) and time (l), Maxwell's equations become
$$ \\begin{align*}
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
\\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &= -\\tilde{\\partial}_t \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}
Poynting vector doc updates
2020-01-12 22:50:01 -08:00
- \\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}} \\\\
lots more docs
2019-12-01 02:32:31 -08:00
\\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*} $$
with
$$ \\begin{align*}
\\hat{B}_\\vec{r} &= \\mu_{\\vec{r} + \\frac{1}{2}} \\cdot \\hat{H}_{\\vec{r} + \\frac{1}{2}} \\\\
\\tilde{D}_\\vec{r} &= \\epsilon_\\vec{r} \\cdot \\tilde{E}_\\vec{r}
\\end{align*} $$
where the spatial subscripts are abbreviated as \\( \\vec{r} = (m, n, p) \\) and
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
\\( \\vec{r} + \\frac{1}{2} = (m + \\frac{1}{2}, n + \\frac{1}{2}, p + \\frac{1}{2}) \\),
\\( \\tilde{E} \\) and \\( \\hat{H} \\) are the electric and magnetic fields,
\\( \\tilde{J} \\) and \\( \\hat{M} \\) are the electric and magnetic current distributions,
and \\( \\epsilon \\) and \\( \\mu \\) are the dielectric permittivity and magnetic permeability.
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:
Poynting vector doc updates
2020-01-12 22:50:01 -08:00
[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]
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
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:
Put H at vertices and label them +1/2
2019-12-10 01:14:21 -08:00
[figure: Field components]
(m - 1/2,=> ____________Hx__________[H] <= r + 1/2 = (m + 1/2,
n + 1/2, /: /: /| n + 1/2,
z y p + 1/2) / : / : / | p + 1/2)
|/_x / : / : / |
/ : Ez__________Hy | Locations of the E- and
/ : : : /| | H-field components for the
(m - 1/2, / : : Ey...../.|..Hz [E] fore-vector at r = (m,n,p)
n - 1/2, =>/________________________/ | /| (the large cube's center)
p + 1/2) | : : / | | / | and [H] back-vector at r + 1/2
| : :/ | |/ | (the top right corner)
| : [E].......|.Ex |
| :.................|......| <= (m + 1/2, n + 1/2, p + 1/2)
| / | /
| / | /
| / | / This is the Yee discretization
| / | / scheme ("Yee cell").
r - 1/2 = | / | /
(m - 1/2, |/ |/
n - 1/2,=> |________________________| <= (m + 1/2, n - 1/2, p - 1/2)
p - 1/2)
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
Each component forms its own grid, offset from the others:
[figure: E-fields for adjacent cells]
Put H at vertices and label them +1/2
2019-12-10 01:14:21 -08:00
H1__________Hx0_________H0
z y /: /|
|/_x / : / | This figure shows H back-vector locations
/ : / | H0, H1, etc. and their associated components
Hy1 : Hy0 | H0 = (Hx0, Hy0, Hz0) etc.
/ : / |
/ Hz1 / Hz0
H2___________Hx3_________H3 | The equivalent drawing for E would have
| : | | fore-vectors located at the cube's
| : | | center (and the centers of adjacent cubes),
| : | | with components on the cube's faces.
| H5..........Hx4...|......H4
| / | /
Hz2 / Hz2 /
| / | /
| Hy6 | Hy4
| / | /
|/ |/
H6__________Hx7__________H7
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
lots more docs
2019-12-01 02:32:31 -08:00
The divergence equations can be derived by taking the divergence of the curl equations
and combining them with charge continuity,
$$ \\hat{\\nabla} \\cdot \\tilde{J} + \\hat{\\partial}_t \\rho = 0 $$
implying that the discrete Maxwell's equations do not produce spurious charges.
Wave equation
-------------
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
Taking the backward curl of the \\( \\tilde{\\nabla} \\times \\tilde{E} \\) equation and
replacing the resulting \\( \\hat{\\nabla} \\times \\hat{H} \\) term using its respective equation,
and setting \\( \\hat{M} \\) to zero, we can form the discrete wave equation:
lots more docs
2019-12-01 02:32:31 -08:00
$$
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
\\begin{align*}
\\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &=
-\\tilde{\\partial}_t \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}
Poynting vector doc updates
2020-01-12 22:50:01 -08:00
- \\hat{M}_{l-1, \\vec{r} + \\frac{1}{2}} \\\\
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
\\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}}) &=
\\hat{\\nabla} \\times (-\\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}}) &=
-\\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}} \\\\
more documentation
2020-01-08 00:51:56 -08:00
\\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}})
lots more docs
2019-12-01 02:32:31 -08:00
+ \\tilde{\\partial}_t \\hat{\\partial}_t \\epsilon_\\vec{r} \\cdot \\tilde{E}_{l, \\vec{r}}
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
&= \\tilde{\\partial}_t \\tilde{J}_{l - \\frac{1}{2}, \\vec{r}}
\\end{align*}
$$
lots more docs
2019-12-01 02:32:31 -08:00
more documentation
2020-01-08 00:51:56 -08:00
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
$$
Poynting vector doc updates
2020-01-12 22:50:01 -08:00
where \\( c = \\sqrt{\\mu \\epsilon} \\).
more documentation
2020-01-08 00:51:56 -08:00
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}) $$
Poynting vector doc updates
2020-01-12 22:50:01 -08:00
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).
more documentation
2020-01-08 00:51:56 -08:00
lots more docs
2019-12-01 02:32:31 -08:00
Grid description
================
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
work on grid text
2019-12-10 01:52:07 -08:00
As described in the section on scalar discrete derivatives above, cell widths
(`dx[i]`, `dy[j]`, `dz[k]`) along each axis can be arbitrary and independently
defined. Moreover, all field components are actually defined at "derived" or "dual"
positions, in-between the "base" grid points on one or more axes.
To get a better sense of how this works, let's start by drawing a grid with uniform
`dy` and `dz` and nonuniform `dx`. We will only draw one cell in the y and z dimensions
to make the illustration simpler; we need at least two cells in the x dimension to
demonstrate how nonuniform `dx` affects the various components.
Place the E fore-vectors at integer indices \\( r = (m, n, p) \\) and the H back-vectors
grid figure fixup
2020-01-04 18:19:04 -08:00
at fractional indices \\( r + \\frac{1}{2} = (m + \\frac{1}{2}, n + \\frac{1}{2},
p + \\frac{1}{2}) \\). Remember that these are indices and not coordinates; they can
correspond to arbitrary (monotonically increasing) coordinates depending on the cell widths.
work on grid text
2019-12-10 01:52:07 -08:00
Draw lines to denote the planes on which the H components and back-vectors are defined.
For simplicity, don't draw the equivalent planes for the E components and fore-vectors,
except as necessary to show their locations -- it's easiest to just connect them to their
grid figure fixup
2020-01-04 18:19:04 -08:00
associated H-equivalents.
The result looks something like this:
more grid work
2019-12-09 21:28:26 -08:00
[figure: Component centers]
grid figure fixup
2020-01-04 18:19:04 -08:00
p=
[H]__________Hx___________[H]_____Hx______[H] __ +1/2
z y /: /: /: /: /| | |
|/_x / : / : / : / : / | | |
/ : / : / : / : / | | |
Hy : Ez...........Hy : Ez......Hy | | |
/: : : : /: : : : /| | | |
/ : Hz : Ey....../.:..Hz : Ey./.|..Hz __ 0 | dz[0]
/ : /: : / / : /: : / / | /| | |
/_________________________/_______________/ | / | | |
| :/ : :/ | :/ : :/ | |/ | | |
| Ex : [E].......|..Ex : [E]..|..Ex | | |
| : | : | | | |
| [H]..........Hx....|......[H].....H|x.....[H] __ --------- (n=+1/2, p=-1/2)
| / | / | / / /
Hz / Hz / Hz / / /
| / | / | / / /
| Hy | Hy | Hy __ 0 / dy[0]
| / | / | / / /
| / | / | / / /
|/ |/ |/ / /
[H]__________Hx___________[H]_____Hx______[H] __ -1/2 /
work on grid text
2019-12-10 01:52:07 -08:00
=n
grid figure fixup
2020-01-04 18:19:04 -08:00
|------------|------------|-------|-------|
-1/2 0 +1/2 +1 +3/2 = m
work on grid text
2019-12-10 01:52:07 -08:00
------------------------- ----------------
dx[0] dx[1]
more grid work
2019-12-09 21:28:26 -08:00
Part of a nonuniform "base grid", with labels specifying
Put H at vertices and label them +1/2
2019-12-10 01:14:21 -08:00
positions of the various field components. [E] fore-vectors
are at the cell centers, and [H] back-vectors are at the
vertices. H components along the near (-y) top (+z) edge
have been omitted to make the insides of the cubes easier
to visualize.
more explanation about grid figures
2020-01-04 18:19:49 -08:00
The above figure shows where all the components are located; however, it is also useful to show
what volumes those components correspond to. Consider the Ex component at `m = +1/2`: it is
shifted in the x-direction by a half-cell from the E fore-vector at `m = 0` (labeled `[E]`
in the figure). It corresponds to a volume between `m = 0` and `m = +1` (the other
dimensions are not shifted, i.e. they are still bounded by `n, p = +-1/2`). (See figure
below). Since `m` is an index and not an x-coordinate, the Ex component is not necessarily
at the center of the volume it represents, and the x-length of its volume is the derived
quantity `dx'[0] = (dx[0] + dx[1]) / 2` rather than the base `dx`.
(See also `Scalar derivatives and cell shifts`).
work on grid text
2019-12-10 01:52:07 -08:00
[figure: Ex volumes]
grid figure fixup
2020-01-04 18:19:04 -08:00
p=
<_________________________________________> __ +1/2
z y << /: / /: >> | |
|/_x < < / : / / : > > | |
< < / : / / : > > | |
< < / : / / : > > | |
<: < / : : / : >: > | |
< : < / : : / : > : > __ 0 | dz[0]
< : < / : : / :> : > | |
<____________/____________________/_______> : > | |
< : < | : : | > : > | |
< Ex < | : Ex | > Ex > | |
< : < | : : | > : > | |
< : <....|.......:........:...|.......>...:...> __ --------- (n=+1/2, p=-1/2)
< : < | / : /| /> : > / /
< : < | / : / | / > : > / /
< :< | / :/ | / > :> / /
< < | / : | / > > _ 0 / dy[0]
< < | / | / > > / /
< < | / | / > > / /
<< |/ |/ >> / /
<____________|____________________|_______> __ -1/2 /
=n
|------------|------------|-------|-------|
-1/2 0 +1/2 +1 +3/2 = m
~------------ -------------------- -------~
dx'[-1] dx'[0] dx'[1]
Put H at vertices and label them +1/2
2019-12-10 01:14:21 -08:00
The Ex values are positioned on the x-faces of the base
grid. They represent the Ex field in volumes shifted by
a half-cell in the x-dimension, as shown here. Only the
grid figure fixup
2020-01-04 18:19:04 -08:00
center cell (with width dx'[0]) is fully shown; the
other two are truncated (shown using >< markers).
Put H at vertices and label them +1/2
2019-12-10 01:14:21 -08:00
Note that the Ex positions are the in the same positions
as the previous figure; only the cell boundaries have moved.
Also note that the points at which Ex is defined are not
necessarily centered in the volumes they represent; non-
uniform cell sizes result in off-center volumes like the
center cell here.
more explanation about grid figures
2020-01-04 18:19:49 -08:00
The next figure shows the volumes corresponding to the Hy components, which
are shifted in two dimensions (x and z) compared to the base grid.
work on grid text
2019-12-10 01:52:07 -08:00
[figure: Hy volumes]
grid figure fixup
2020-01-04 18:19:04 -08:00
p=
z y mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm __ +1/2 s
|/_x << m: m: >> | |
< < m : m : > > | | dz'[1]
< < m : m : > > | |
Hy........... m........Hy...........m......Hy > | |
< < m : m : > > | |
< ______ m_____:_______________m_____:_>______ __ 0
< < m /: m / > > | |
mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm > | |
< < | / : | / > > | | dz'[0]
< < | / : | / > > | |
< < | / : | / > > | |
< wwwww|w/wwwwwwwwwwwwwwwwwww|w/wwwww>wwwwwwww __ s
< < |/ w |/ w> > / /
_____________|_____________________|________ > / /
< < | w | w > > / /
< Hy........|...w........Hy.......|...w...>..Hy _ 0 / dy[0]
< < | w | w > > / /
<< | w | w > > / /
< |w |w >> / /
wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww __ -1/2 /
|------------|------------|--------|-------|
-1/2 0 +1/2 +1 +3/2 = m
Put H at vertices and label them +1/2
2019-12-10 01:14:21 -08:00
~------------ --------------------- -------~
grid figure fixup
2020-01-04 18:19:04 -08:00
dx'[-1] dx'[0] dx'[1]
Put H at vertices and label them +1/2
2019-12-10 01:14:21 -08:00
The Hy values are positioned on the y-edges of the base
grid. Again here, the 'Hy' labels represent the same points
as in the basic grid figure above; the edges have shifted
by a half-cell along the x- and z-axes.
The grid lines _|:/ are edges of the area represented by
each Hy value, and the lines drawn using <m>.w represent
edges where a cell's faces extend beyond the drawn area
(i.e. where the drawing is truncated in the x- or z-
directions).
more grid work
2019-12-09 21:28:26 -08:00
lots more fdmath documentation
2019-12-08 01:46:47 -08:00
explain dx_lists_t
2020-01-04 18:19:57 -08:00
Datastructure: dx_lists_t
-------------------
move stuff under fdmath
2019-11-27 22:59:52 -08:00
explain dx_lists_t
2020-01-04 18:19:57 -08:00
In this documentation, the E fore-vectors are placed on the base grid. An
equivalent formulation could place the H back-vectors on the base grid instead.
However, in the case of a non-uniform grid, the operation to get from the "base"
cell widths to "derived" ones is not its own inverse.
work on grid text
2019-12-10 01:52:07 -08:00
explain dx_lists_t
2020-01-04 18:19:57 -08:00
The base grid's cell sizes could be fully described by a list of three 1D arrays,
specifying the cell widths along all three axes:
[dx, dy, dz] = [[dx[0], dx[1], ...], [dy[0], ...], [dz[0], ...]]
Note that this is a list-of-arrays rather than a 2D array, as the simulation domain
may have a different number of cells along each axis.
Knowing the base grid's cell widths and the boundary conditions (periodic unless
otherwise noted) is enough information to calculate the cell widths `dx'`, `dy'`,
and `dz'` for the derived grids.
However, since most operations are trivially generalized to allow either E or H
to be defined on the base grid, they are written to take the a full set of base
and derived cell widths, distinguished by which field they apply to rather than
their "base" or "derived" status. This removes the need for each function to
generate the derived widths, and makes the "base" vs "derived" distinction
unnecessary in the code.
The resulting data structure containing all the cell widths takes the form of a
list-of-lists-of-arrays. The first list-of-arrays provides the cell widths for
the E-field fore-vectors, while the second list-of-arrays does the same for the
H-field back-vectors:
[[[dx_e[0], dx_e[1], ...], [dy_e[0], ...], [dz_e[0], ...]],
[[dx_h[0], dx_h[1], ...], [dy_h[0], ...], [dz_h[0], ...]]]
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.
work on grid text
2019-12-10 01:52:07 -08:00
more documentation
2020-01-08 00:51:56 -08:00
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.
documentation updates
2019-11-26 01:47:52 -08:00
"""
move stuff under fdmath
2019-11-27 22:59:52 -08:00
from .types import fdfield_t, vfdfield_t, dx_lists_t, fdfield_updater_t
from .vectorization import vec, unvec
from . import operators, functional, types, vectorization
Reference in New Issue Copy Permalink