diff --git a/meanas/fdfd/__init__.py b/meanas/fdfd/__init__.py index 240d711..e94adc3 100644 --- a/meanas/fdfd/__init__.py +++ b/meanas/fdfd/__init__.py @@ -19,39 +19,38 @@ Submodules: ================================================================ -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: +From the "Frequency domain" section of `meanas.fdmath`, we have $$ \begin{aligned} - \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} \\ + \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} \\ \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t \end{aligned} $$ -The difference multipliers for an integer-indexed sequence are +resulting in $$ \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} $$ -The half-step factors cancel when the staggered sample times are included. -Maxwell's equations for the phasors are then +Maxwell's equations are then $$ \begin{aligned} \tilde{\nabla} \times \tilde{E}_{\vec{r}} &= - -\imath \Omega \hat{B}_{\vec{r} + \frac{1}{2}} + \imath \Omega e^{-\imath \omega \Delta_t / 2} \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 e^{ \imath \omega \Delta_t / 2} \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}} @@ -62,9 +61,13 @@ 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} $$ @@ -73,17 +76,14 @@ 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 ce9827b..857cf18 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^x_{m + \frac{1}{2},n,p}) + &+ \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{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, \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}}) &= @@ -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}} - \tilde{\partial}_t \tilde{J}_{l-\frac{1}{2},\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}}) + \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,42 +403,37 @@ $$ Frequency domain ---------------- -We use the $e^{+\imath \omega t}$ convention, matching the FDFD operators. -Substitute time-harmonic fields with phases referenced to their actual sample times: +We can substitute in a time-harmonic fields $$ \begin{aligned} - \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} + \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} \end{aligned} $$ -For an integer-indexed sequence proportional to $e^{+\imath \omega l \Delta_t}$, -the forward and backward difference multipliers are +resulting in $$ \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} $$ -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: +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}} \\ + -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\ $$ Plane waves and Dispersion relation ------------------------------------ -With uniform grid spacings, a uniform material distribution, and no sources +With uniform material distribution and no sources $$ \begin{aligned} @@ -479,31 +474,31 @@ resulting in $$ \begin{aligned} - \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}\\ + \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{aligned} $$ -with similar expressions for the y and z dimensions (and $K_y, K_z$). +with similar expressions for the y and z dimnsions (and $K_y, K_z$). This implies $$ - \tilde{\nabla}^2 \phi_{\vec{r}} = -(K_x^2 + K_y^2 + K_z^2) \phi_{\vec{r}} \\ + \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 $$ -where $c = 1 / \sqrt{\mu \epsilon}$. +where $c = \sqrt{\mu \epsilon}$. -To keep $\omega$ real for every real grid wavevector $(k_x, k_y, k_z)$, -we use the conservative stability bound +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 < \Delta_x / \sqrt{3}$. -This limits the distance that light travels in one timestep to less than -$\Delta_x / \sqrt{3}$ on a uniform cubic grid. +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/fdmath/functional.py b/meanas/fdmath/functional.py index c48848a..dacde37 100644 --- a/meanas/fdmath/functional.py +++ b/meanas/fdmath/functional.py @@ -3,6 +3,7 @@ Math functions for finite difference simulations Basic discrete calculus etc. """ +from typing import TypeVar from collections.abc import Sequence, Callable import numpy @@ -31,15 +32,14 @@ def _derivative(axis: int, forward: bool, widths: NDArray | None) -> fdfield_upd raise ValueError('out must not alias the field or widths') if not output.flags.writeable or not numpy.can_cast(dtype, output.dtype, casting='same_kind'): raise ValueError('out must be writable with a compatible dtype') - dtype = numpy.result_type(dtype, output.dtype) low, high = [slice(None)] * 3, [slice(None)] * 3 low[axis], high[axis] = slice(None, -1), slice(1, None) first, last = [slice(None)] * 3, [slice(None)] * 3 first[axis], last[axis] = slice(0, 1), slice(-1, None) destination = low if forward else high - numpy.subtract(f[tuple(high)], f[tuple(low)], out=output[tuple(destination)], dtype=dtype) + numpy.subtract(f[tuple(high)], f[tuple(low)], out=output[tuple(destination)]) destination = last if forward else first - numpy.subtract(f[tuple(first)], f[tuple(last)], out=output[tuple(destination)], dtype=dtype) + numpy.subtract(f[tuple(first)], f[tuple(last)], out=output[tuple(destination)]) if widths is not None: expand = [1, 1, 1] expand[axis] = len(widths) @@ -55,7 +55,6 @@ def deriv_forward( Each callable accepts ``derivative(field, *, out=None)``. An output buffer must match the field shape and must not alias the field or width arrays. - Subtraction uses the dtype promoted by the field, widths, and output buffer. Slice-based subtraction avoids a full-size roll temporary. """ return tuple(_derivative(axis, True, None if dx_e is None else dx_e[axis]) for axis in range(3)) @@ -71,9 +70,12 @@ def deriv_back( return tuple(_derivative(axis, False, None if dx_h is None else dx_h[axis]) for axis in range(3)) +TT = TypeVar('TT', bound='NDArray[floating | complexfloating]') + + def curl_forward( dx_e: Sequence[NDArray[floating | complexfloating]] | None = None, - ) -> Callable[[NDArray[floating | complexfloating]], NDArray[floating | complexfloating]]: + ) -> Callable[[TT], TT]: r""" Curl operator for use with the E field. @@ -83,14 +85,12 @@ def curl_forward( Returns: Function `f` for taking the discrete forward curl of a field, - `f(E)` -> curlE $= \nabla_f \times E$. - The output dtype is promoted from the field and cell widths. + `f(E)` -> curlE $= \nabla_f \times E$ """ Dx, Dy, Dz = deriv_forward(dx_e) - def ce_fun(e: NDArray[floating | complexfloating]) -> NDArray[floating | complexfloating]: - dtype = e.dtype if dx_e is None else numpy.result_type(e, *dx_e) - output = numpy.empty_like(e, dtype=dtype) + def ce_fun(e: TT) -> TT: + output = numpy.empty_like(e) output[0] = Dy(e[2]) output[1] = Dz(e[0]) output[2] = Dx(e[1]) @@ -104,7 +104,7 @@ def curl_forward( def curl_back( dx_h: Sequence[NDArray[floating | complexfloating]] | None = None, - ) -> Callable[[NDArray[floating | complexfloating]], NDArray[floating | complexfloating]]: + ) -> Callable[[TT], TT]: r""" Create a function which takes the backward curl of a field. @@ -114,14 +114,12 @@ def curl_back( Returns: Function `f` for taking the discrete backward curl of a field, - `f(H)` -> curlH $= \nabla_b \times H$. - The output dtype is promoted from the field and cell widths. + `f(H)` -> curlH $= \nabla_b \times H$ """ Dx, Dy, Dz = deriv_back(dx_h) - def ch_fun(h: NDArray[floating | complexfloating]) -> NDArray[floating | complexfloating]: - dtype = h.dtype if dx_h is None else numpy.result_type(h, *dx_h) - output = numpy.empty_like(h, dtype=dtype) + def ch_fun(h: TT) -> TT: + output = numpy.empty_like(h) output[0] = Dy(h[2]) output[1] = Dz(h[0]) output[2] = Dx(h[1]) diff --git a/meanas/fdmath/operators.py b/meanas/fdmath/operators.py index d9348ba..8b7cabc 100644 --- a/meanas/fdmath/operators.py +++ b/meanas/fdmath/operators.py @@ -61,9 +61,6 @@ def shift_with_mirror( Utility operator for performing an n-element shift along a specified axis, with mirror boundary conditions applied to the cells beyond the receding edge. - The absolute shift distance may not exceed the axis length. A shift by the - full axis length reverses that axis; shifts on a singleton axis are identities. - Args: axis: Axis to shift along. x=0, y=1, z=2 shape: Shape of the grid being shifted @@ -76,8 +73,8 @@ def shift_with_mirror( raise Exception(f'Invalid shape: {shape}') if axis not in range(len(shape)): raise Exception(f'Invalid direction: {axis}, shape is {shape}') - if abs(shift_distance) > shape[axis]: - raise ValueError(f'Shift ({shift_distance}) is too large for axis {axis} of size {shape[axis]}') + if shift_distance >= shape[axis]: + raise Exception(f'Shift ({shift_distance}) is too large for axis {axis} of size {shape[axis]}') def mirrored_range(n: int, s: int) -> NDArray[numpy.int_]: v = numpy.arange(n) + s diff --git a/meanas/fdmath/vectorization.py b/meanas/fdmath/vectorization.py index 2a47fd8..77b722f 100644 --- a/meanas/fdmath/vectorization.py +++ b/meanas/fdmath/vectorization.py @@ -113,7 +113,7 @@ def unvec( Returns `None` if called with `v=None`. Args: - v: 1D array-like representing a vector field of shape shape (or None) + v: 1D ndarray representing a vector field of shape shape (or None) shape: shape of the vector field nvdim: Number of components in each vector @@ -122,4 +122,4 @@ def unvec( """ if v is None: return None - return numpy.asarray(v).reshape((nvdim, *shape), order='C') + return v.reshape((nvdim, *shape), order='C') # type: ignore diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py index 769f362..480ed87 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 += dt * curl_back(H[t]) / epsilon + E += curl_back(H[t]) / epsilon The full update should be - E += dt * (curl_back(H[t]) - J) / epsilon + E += (curl_back(H[t]) + J) / epsilon - which requires an additional step of `E -= dt * J / epsilon` which is not performed + which requires an additional step of `E += 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 -= dt * curl_forward(E[t]) / mu + H -= curl_forward(E[t]) / mu The full update should be - H -= dt * (curl_forward(E[t]) + M) / mu + H -= (curl_forward(E[t]) + M) / mu - which requires an additional step of `H -= dt * 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`. diff --git a/meanas/test/test_fdmath_functional.py b/meanas/test/test_fdmath_functional.py index 3a1f3cf..f6d1491 100644 --- a/meanas/test/test_fdmath_functional.py +++ b/meanas/test/test_fdmath_functional.py @@ -62,55 +62,6 @@ def test_curls_match_sparse_operators_on_nonuniform_grid() -> None: assert_fields_close(curl_back, matrix_back, atol=1e-12, rtol=1e-12) -@pytest.mark.parametrize('name', ['curl_forward', 'curl_back']) -@pytest.mark.parametrize('field_dtype', [numpy.float32, numpy.complex64]) -@pytest.mark.parametrize('complex_widths', [False, True]) -def test_curls_promote_field_dtype_and_match_sparse_operators(name, field_dtype, complex_widths) -> None: - field = VECTOR_FIELD.astype(field_dtype) if numpy.issubdtype(field_dtype, numpy.complexfloating) else VECTOR_FIELD.real.astype(field_dtype) - widths = [dx * (1 + 0.5j) if complex_widths else dx for dx in DX_E] - result = getattr(fd_functional, name)(widths)(field) - expected = unvec(getattr(fd_operators, name)(widths) @ vec(field), SHAPE) - assert result.dtype == numpy.result_type(field, *widths) - assert_fields_close(result, expected, atol=1e-12, rtol=1e-12) - - -@pytest.mark.parametrize('name', ['curl_forward', 'curl_back']) -def test_curls_preserve_small_residual_between_large_derivative_terms(name) -> None: - shape = (1, 2, 2) - field = numpy.zeros((3, *shape), dtype=numpy.float32) - field[2, :, 1, :] = 1e8 - field[1, :, :, 1] = 1e8 - widths = [numpy.ones(1), numpy.full(2, 1e8 / (1e8 + 1)), numpy.ones(2)] - result = getattr(fd_functional, name)(widths)(field) - expected = unvec(getattr(fd_operators, name)(widths) @ vec(field), shape) - assert abs(expected[0, 0, 0, 0]) == pytest.approx(1) - assert_fields_close(result, expected, atol=1e-7, rtol=1e-12) - - -@pytest.mark.parametrize('name', ['deriv_forward', 'deriv_back']) -@pytest.mark.parametrize('axis', [0, 1, 2]) -@pytest.mark.parametrize('spaced', [False, True]) -def test_derivatives_promote_subtraction_before_overflow(name, axis, spaced) -> None: - shape = [1, 1, 1] - shape[axis] = 2 - field = numpy.array([3e38, -3e38], dtype=numpy.float32).reshape(shape) - widths = [numpy.full(n, 1e38) for n in shape] if spaced else None - derivative = getattr(fd_functional, name)(widths)[axis] - wide_field = field.astype(numpy.float64) - if name == 'deriv_forward': - expected = numpy.roll(wide_field, -1, axis=axis) - wide_field - else: - expected = wide_field - numpy.roll(wide_field, 1, axis=axis) - if spaced: - expected /= 1e38 - with numpy.errstate(over='raise', invalid='raise'): - if spaced: - assert_close(derivative(field), expected) - output = numpy.empty(shape, dtype=numpy.float64) - assert derivative(field, out=output) is output - assert_close(output, expected) - - @pytest.mark.parametrize('periodic', [True, False]) def test_gridlock_spacings_differentiate_physical_coordinates(periodic) -> None: grid = Grid([[0, 1, 3, 7, 15], [-2, 0, 5, 6], [3, 6, 10]], periodic=periodic) diff --git a/meanas/test/test_fdmath_operators.py b/meanas/test/test_fdmath_operators.py index 2ba7755..9b1dec0 100644 --- a/meanas/test/test_fdmath_operators.py +++ b/meanas/test/test_fdmath_operators.py @@ -56,8 +56,7 @@ def test_shift_circ_rejects_invalid_arguments(args: tuple[int, tuple[int, ...], [ ((0, (2,), 1), 'Invalid shape'), ((3, SHAPE, 1), 'Invalid direction'), - ((0, SHAPE, SHAPE[0] + 1), 'too large'), - ((0, SHAPE, -SHAPE[0] - 1), 'too large'), + ((0, SHAPE, SHAPE[0]), 'too large'), ], ) def test_shift_with_mirror_rejects_invalid_arguments(args: tuple[int, tuple[int, ...], int], message: str) -> None: @@ -65,23 +64,6 @@ def test_shift_with_mirror_rejects_invalid_arguments(args: tuple[int, tuple[int, operators.shift_with_mirror(*args) -@pytest.mark.parametrize('axis', [0, 1, 2]) -@pytest.mark.parametrize('polarity', [-1, 1]) -def test_shift_with_mirror_by_full_axis_length_reverses_axis(axis, polarity) -> None: - result = _apply_scalar_matrix(operators.shift_with_mirror(axis, SHAPE, polarity * SHAPE[axis])) - assert_close(result, numpy.flip(SCALAR_FIELD, axis=axis)) - - -@pytest.mark.parametrize('axis', [0, 1, 2]) -@pytest.mark.parametrize('shift_distance', [-1, 0, 1]) -def test_shift_with_mirror_on_singleton_axis_is_identity(axis, shift_distance) -> None: - shape = list(SHAPE) - shape[axis] = 1 - field = numpy.arange(numpy.prod(shape)).reshape(shape) - result = (operators.shift_with_mirror(axis, shape, shift_distance) @ field.ravel()).reshape(shape) - assert_close(result, field) - - def test_vec_cross_matches_pointwise_cross_product() -> None: matrix_result = unvec(operators.vec_cross(vec(VECTOR_LEFT)) @ vec(VECTOR_RIGHT), SHAPE) expected = numpy.empty_like(VECTOR_LEFT) diff --git a/meanas/test/test_fdmath_vectorization.py b/meanas/test/test_fdmath_vectorization.py index a895d24..c55f7cd 100644 --- a/meanas/test/test_fdmath_vectorization.py +++ b/meanas/test/test_fdmath_vectorization.py @@ -1,5 +1,4 @@ import numpy -import pytest from ..fdmath import unvec, vec from ._test_builders import complex_ramp, real_ramp @@ -45,11 +44,3 @@ def test_unvec_with_two_components_round_trips_vector() -> None: def test_vec_accepts_arraylike_input() -> None: arraylike = [[[1, 2], [3, 4]], [[5, 6], [7, 8]]] assert_close(vec(arraylike), numpy.ravel(arraylike, order='C')) - - -@pytest.mark.parametrize('arraylike_type', [list, tuple]) -def test_unvec_accepts_arraylike_input(arraylike_type) -> None: - vector = numpy.arange(3 * numpy.prod(SHAPE)) * (1 + 0.5j) - field = unvec(arraylike_type(vector), SHAPE) - assert field.shape == (3, *SHAPE) - assert_close(vec(field), vector) diff --git a/meanas/test/test_numerical_feedback.py b/meanas/test/test_numerical_feedback.py index 6c2d27b..1658bb7 100644 --- a/meanas/test/test_numerical_feedback.py +++ b/meanas/test/test_numerical_feedback.py @@ -289,10 +289,7 @@ def test_slice_derivatives_match_periodic_reference_and_output_buffers(dtype, ax f = np.ascontiguousarray(f) widths = [np.linspace(.7, 1.3, n) for n in f.shape] if spaced else None derivative = (deriv_forward if forward else deriv_back)(widths)[axis] - # Widths promote the calculation before subtraction, including integer fields. - calculation_field = f.astype(np.result_type(f, widths[axis])) if spaced else f - expected = (np.roll(calculation_field, -1, axis=axis) - calculation_field if forward - else calculation_field - np.roll(calculation_field, 1, axis=axis)) + expected = np.roll(f, -1, axis=axis) - f if forward else f - np.roll(f, 1, axis=axis) if spaced: shape = [1, 1, 1] shape[axis] = f.shape[axis]