diff --git a/README.md b/README.md index 3dcdfe0..c3632c0 100644 --- a/README.md +++ b/README.md @@ -94,6 +94,59 @@ python3 -m pytest -rsxX | tee test_results.txt ## Use -See `examples/` for some simple examples; you may need additional -packages such as [gridlock](https://mpxd.net/code/jan/gridlock) -to run the examples. +`meanas` is organized around a few core workflows: + +- `meanas.fdfd`: frequency-domain wave equations, sparse operators, SCPML, and + iterative solves for driven problems. +- `meanas.fdfd.waveguide_2d` / `meanas.fdfd.waveguide_3d`: waveguide mode + solvers, mode-source construction, and overlap windows for port-based + excitation and analysis. +- `meanas.fdtd`: Yee-step updates, CPML boundaries, flux/energy accounting, and + on-the-fly phasor extraction for comparing time-domain runs against FDFD. +- `meanas.fdmath`: low-level finite-difference operators, vectorization helpers, + and derivations shared by the FDTD and FDFD layers. + +The most mature user-facing workflows are: + +1. Build an FDFD operator or waveguide port source, then solve a driven + frequency-domain problem. +2. Run an FDTD simulation, extract one or more frequency-domain phasors with + `meanas.fdtd.accumulate_phasor(...)`, and compare those phasors against an + FDFD reference on the same Yee grid. + +Tracked examples under `examples/` are the intended starting points: + +- `examples/fdtd.py`: broadband FDTD pulse excitation, phasor extraction, and a + residual check against the matching FDFD operator. +- `examples/waveguide.py`: waveguide mode solving, unidirectional mode-source + construction, overlap readout, and FDTD/FDFD comparison on a guided structure. +- `examples/fdfd.py`: direct frequency-domain waveguide excitation and overlap / + Poynting analysis without a time-domain run. + +Several examples rely on optional packages such as +[gridlock](https://mpxd.net/code/jan/gridlock). + +### Frequency-domain waveguide workflow + +For a structure with a constant cross-section in one direction: + +1. Build `dxes` and the diagonal `epsilon` / `mu` distributions on the Yee grid. +2. Solve the port mode with `meanas.fdfd.waveguide_3d.solve_mode(...)`. +3. Build a unidirectional source with `compute_source(...)`. +4. Build a matching overlap window with `compute_overlap_e(...)`. +5. Solve the full FDFD problem and project the result onto the overlap window or + evaluate plane flux with `meanas.fdfd.functional.poynting_e_cross_h(...)`. + +### Time-domain phasor workflow + +For a broadband or continuous-wave FDTD run: + +1. Advance the fields with `meanas.fdtd.maxwell_e/maxwell_h` or + `updates_with_cpml(...)`. +2. Inject electric current using the same sign convention used throughout the + examples and library: `E -= dt * J / epsilon`. +3. Accumulate the desired phasor with `accumulate_phasor(...)` or the Yee-aware + wrappers `accumulate_phasor_e/h/j(...)`. +4. Build the matching FDFD operator on the stretched `dxes` if CPML/SCPML is + part of the simulation, and compare the extracted phasor to the FDFD field or + residual. diff --git a/examples/fdtd.py b/examples/fdtd.py index 2a50ddc..704c2f1 100644 --- a/examples/fdtd.py +++ b/examples/fdtd.py @@ -1,7 +1,12 @@ """ -Example code for running an FDTD simulation +Example code for a broadband FDTD run with phasor extraction. -See main() for simulation setup. +This script shows the intended low-level workflow for: + +1. building a Yee-grid simulation with CPML on all faces, +2. driving it with an electric-current pulse, +3. extracting a single-frequency phasor on the fly, and +4. checking that phasor against the matching stretched-grid FDFD operator. """ import sys @@ -150,7 +155,8 @@ def main(): # print(f'Save time interval would be {sample_interval} * dt = {sample_interval * dt:3g}') - # Source parameters and function + # Source parameters and function. The pulse normalization is kept outside + # accumulate_phasor(); the helper only performs the Fourier sum. source_phasor, delay = gaussian_packet(wl=wl, dwl=100, dt=dt, turn_on=1e-5) aa, cc, ss = source_phasor(numpy.arange(max_t)) srca_real = aa * cc @@ -170,7 +176,8 @@ def main(): update_E(ee, hh, epsilon) if tt < src_maxt: - # This codebase uses E -= dt * J / epsilon for electric-current injection. + # Electric-current injection uses E -= dt * J / epsilon, which is + # the same sign convention used by the matching FDFD right-hand side. ee[1, *(grid.shape // 2)] -= srca_real[tt] update_H(ee, hh) @@ -193,9 +200,11 @@ def main(): dt, ee, tt, - # The pulse is delayed relative to t=0, so the readout needs the same phase shift. + # The pulse is delayed relative to t=0, so the extracted phasor + # needs the same phase offset in its sample times. offset_steps=0.5 - delay / dt, - # accumulate_phasor() already includes dt, so undo the dt in phasor_norm here. + # accumulate_phasor() already multiplies by dt, so pass the + # discrete-sum normalization without its extra dt factor. weight=phasor_norm / dt, ) @@ -205,6 +214,8 @@ def main(): for pp in (-1, +1): for dd in range(3): stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_air ** 2, thickness=pml_thickness) + # Compare the extracted phasor to the FDFD operator on the stretched grid, + # not the unstretched Yee spacings used by the raw time-domain update. A = e_full(omega=omega, dxes=dxes_fdfd, epsilon=epsilon) residual = norm(A @ vec(Eph) - vec(b)) / norm(vec(b)) print(f'FDFD residual is {residual}') diff --git a/examples/waveguide.py b/examples/waveguide.py index 1bb02fa..a100ee3 100644 --- a/examples/waveguide.py +++ b/examples/waveguide.py @@ -1,7 +1,11 @@ """ -Example code for running an OpenCL FDTD simulation +Example code for guided-wave FDFD and FDTD comparison. -See main() for simulation setup. +This example is the reference workflow for: + +1. solving a waveguide port mode, +2. turning that mode into a one-sided source and overlap window, +3. comparing a direct FDFD solve against a time-domain phasor extracted from FDTD. """ from typing import Callable import logging @@ -78,7 +82,7 @@ def get_waveguide_mode( omega: float, epsilon: fdfield_t, ) -> tuple[vcfdfield_t, vcfdfield_t]: - """ Create a mode source and overlap window """ + """Create a mode source and overlap window for one forward-going port.""" dims = numpy.array([[-10, -20, -15], [-10, 20, 15]]) * [[numpy.median(numpy.real(dx)) for dx in dxes[0]]] ind_dims = (grid.pos2ind(dims[0], which_shifts=None).astype(int), @@ -94,6 +98,8 @@ def get_waveguide_mode( J = waveguide_3d.compute_source(E=wg_results['E'], wavenumber=wg_results['wavenumber'], omega=omega, epsilon=epsilon, **wg_args) + # compute_overlap_e() returns the normalized upstream overlap window used to + # project another field onto this same guided mode. e_overlap = waveguide_3d.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args, omega=omega) return J, e_overlap @@ -111,7 +117,8 @@ def main( grid, epsilon = draw_grid(dx=dx, pml_thickness=pml_thickness) - # Add PML + # Add SCPML stretching to the FDFD grid; this matches the CPML-backed FDTD + # run below so the two solvers see the same absorbing boundary profile. dxes = [grid.dxyz, grid.autoshifted_dxyz()] for a in (0, 1, 2): for p in (-1, 1): @@ -275,7 +282,8 @@ def main2(): # print(f'Save time interval would be {sample_interval} * dt = {sample_interval * dt:3g}') - # Source parameters and function + # Source parameters and function. The phasor helper only performs the + # Fourier accumulation; the pulse normalization stays explicit here. source_phasor, delay = gaussian_packet(wl=wl, dwl=100, dt=dt, turn_on=1e-5) aa, cc, ss = source_phasor(numpy.arange(max_t)) srca_real = aa * cc @@ -295,7 +303,8 @@ def main2(): update_E(ee, hh, epsilon) if tt < src_maxt: - # This codebase uses E -= dt * J / epsilon for electric-current injection. + # Electric-current injection uses E -= dt * J / epsilon, which is + # the sign convention matched by the FDFD source term -1j * omega * J. ee[1, *(grid.shape // 2)] -= srca_real[tt] update_H(ee, hh) @@ -317,9 +326,11 @@ def main2(): dt, ee, tt, - # The pulse is delayed relative to t=0, so the readout needs the same phase shift. + # The pulse is delayed relative to t=0, so the extracted phasor must + # apply the same delay to its sample times. offset_steps=0.5 - delay / dt, - # accumulate_phasor() already includes dt, so undo the dt in phasor_norm here. + # accumulate_phasor() already contributes dt, so remove the extra dt + # from the externally computed normalization. weight=phasor_norm / dt, ) @@ -329,6 +340,8 @@ def main2(): for pp in (-1, +1): for dd in range(3): stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_cladding ** 2, thickness=pml_thickness) + # Residuals must be checked on the stretched FDFD grid, because the FDTD run + # already includes those same absorbing layers through CPML. A = e_full(omega=omega, dxes=dxes_fdfd, epsilon=epsilon) residual = norm(A @ vec(Eph) - vec(b)) / norm(vec(b)) print(f'FDFD residual is {residual}') diff --git a/meanas/fdfd/__init__.py b/meanas/fdfd/__init__.py index ba57fc4..e94adc3 100644 --- a/meanas/fdfd/__init__.py +++ b/meanas/fdfd/__init__.py @@ -9,9 +9,12 @@ Submodules: - `operators`, `functional`: General FDFD problem setup. - `solvers`: Solver interface and reference implementation. -- `scpml`: Stretched-coordinate perfectly matched layer (scpml) boundary conditions +- `scpml`: Stretched-coordinate perfectly matched layer (SCPML) boundary conditions. - `waveguide_2d`: Operators and mode-solver for waveguides with constant cross-section. -- `waveguide_3d`: Functions for transforming `waveguide_2d` results into 3D. +- `waveguide_3d`: Functions for transforming `waveguide_2d` results into 3D, + including mode-source and overlap-window construction. +- `farfield`, `bloch`, `eme`: specialized helper modules for near/far transforms, + Bloch-periodic problems, and eigenmode expansion. ================================================================ @@ -86,10 +89,6 @@ $$ -\omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \omega \tilde{J}_{\vec{r}} \\ $$ -# TODO FDFD? -# TODO PML - - """ from . import ( solvers as solvers, diff --git a/meanas/fdfd/functional.py b/meanas/fdfd/functional.py index 9d98798..ddd074b 100644 --- a/meanas/fdfd/functional.py +++ b/meanas/fdfd/functional.py @@ -158,10 +158,23 @@ def e_tfsf_source( epsilon: fdfield, mu: fdfield | None = None, ) -> cfdfield_updater_t: - """ + r""" Operator that turns an E-field distribution into a total-field/scattered-field (TFSF) source. + If `A` is the full wave operator from `e_full(...)` and `Q` is the diagonal + mask selecting the total-field region, then the TFSF source is the commutator + + $$ + \frac{A Q - Q A}{-i \omega} E. + $$ + + This vanishes in the interior of the total-field and scattered-field regions + and is supported only at their shared boundary, where the mask discontinuity + makes `A` and `Q` fail to commute. The returned current is therefore the + distributed source needed to inject the desired total field without also + forcing the scattered-field region. + Args: TF_region: mask which is set to 1 in the total-field region, and 0 elsewhere (i.e. in the scattered-field region). @@ -175,7 +188,6 @@ def e_tfsf_source( Function `f` which takes an E field and returns a current distribution, `f(E)` -> `J` """ - # TODO documentation A = e_full(omega, dxes, epsilon, mu) def op(e: cfdfield) -> cfdfield_t: @@ -188,7 +200,13 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield, cfdfield], cfdfi r""" Generates a function that takes the single-frequency `E` and `H` fields and calculates the cross product `E` x `H` = $E \times H$ as required - for the Poynting vector, $S = E \times H$ + for the Poynting vector, $S = E \times H$. + + On the Yee grid, the electric and magnetic components are not stored at the + same locations. This helper therefore applies the same one-cell electric-field + shifts used by the sparse `operators.poynting_e_cross(...)` construction so + that the discrete cross product matches the face-centered energy flux used in + `meanas.fdtd.energy.poynting(...)`. Note: This function also shifts the input `E` field by one cell as required @@ -204,7 +222,8 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield, cfdfield], cfdfi dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` Returns: - Function `f` that returns E x H as required for the poynting vector. + Function `f` that returns the staggered-grid cross product `E \times H`. + For time-average power, call it as `f(E, H.conj())` and take `Re(...) / 2`. """ def exh(e: cfdfield, h: cfdfield) -> cfdfield_t: s = numpy.empty_like(e) diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py index 1282ea6..18aaade 100644 --- a/meanas/fdfd/operators.py +++ b/meanas/fdfd/operators.py @@ -310,16 +310,22 @@ def m2j( def poynting_e_cross(e: vcfdfield, dxes: dx_lists_t) -> sparse.sparray: - """ - Operator for computing the Poynting vector, containing the - (E x) portion of the Poynting vector. + r""" + Operator for computing the staggered-grid `(E \times)` part of the Poynting vector. + + On the Yee grid the E and H components live on different edges, so the + electric field must be shifted by one cell in the appropriate direction + before forming the discrete cross product. This sparse operator encodes that + shifted cross product directly and is the matrix equivalent of + `functional.poynting_e_cross_h(...)`. Args: e: Vectorized E-field for the ExH cross product dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` Returns: - Sparse matrix containing (E x) portion of Poynting cross product. + Sparse matrix containing the `(E \times)` part of the staggered Poynting + cross product. """ shape = [len(dx) for dx in dxes[0]] @@ -339,15 +345,26 @@ def poynting_e_cross(e: vcfdfield, dxes: dx_lists_t) -> sparse.sparray: def poynting_h_cross(h: vcfdfield, dxes: dx_lists_t) -> sparse.sparray: - """ - Operator for computing the Poynting vector, containing the (H x) portion of the Poynting vector. + r""" + Operator for computing the staggered-grid `(H \times)` part of the Poynting vector. + + Together with `poynting_e_cross(...)`, this provides the matrix form of the + Yee-grid cross product used in the flux helpers. The two are related by the + usual antisymmetry of the cross product, + + $$ + H \times E = -(E \times H), + $$ + + once the same staggered field placement is used on both sides. Args: h: Vectorized H-field for the HxE cross product dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` Returns: - Sparse matrix containing (H x) portion of Poynting cross product. + Sparse matrix containing the `(H \times)` part of the staggered Poynting + cross product. """ shape = [len(dx) for dx in dxes[0]] @@ -372,11 +389,23 @@ def e_tfsf_source( epsilon: vfdfield, mu: vfdfield | None = None, ) -> sparse.sparray: - """ + r""" Operator that turns a desired E-field distribution into a total-field/scattered-field (TFSF) source. - TODO: Reference Rumpf paper + Let `A` be the full wave operator from `e_full(...)`, and let + `Q = \mathrm{diag}(TF_region)` be the projector onto the total-field region. + Then the TFSF current operator is the commutator + + $$ + \frac{A Q - Q A}{-i \omega}. + $$ + + Inside regions where `Q` is locally constant, `A` and `Q` commute and the + source vanishes. Only cells at the TF/SF boundary contribute nonzero current, + which is exactly the desired distributed source for injecting the chosen + field into the total-field region without directly forcing the + scattered-field region. Args: TF_region: Mask, which is set to 1 inside the total-field region and 0 in the @@ -388,9 +417,7 @@ def e_tfsf_source( Returns: Sparse matrix that turns an E-field into a current (J) distribution. - """ - # TODO documentation A = e_full(omega, dxes, epsilon, mu) Q = sparse.diags_array(TF_region) return (A @ Q - Q @ A) / (-1j * omega) @@ -404,11 +431,17 @@ def e_boundary_source( mu: vfdfield | None = None, periodic_mask_edges: bool = False, ) -> sparse.sparray: - """ + r""" Operator that turns an E-field distrubtion into a current (J) distribution along the edges (external and internal) of the provided mask. This is just an `e_tfsf_source()` with an additional masking step. + Equivalently, this helper first constructs the TFSF commutator source for the + full mask and then zeroes out all cells except the mask boundary. The + boundary is defined as the set of cells whose mask value changes under a + one-cell shift in any Cartesian direction. With `periodic_mask_edges=False` + the shifts mirror at the domain boundary; with `True` they wrap periodically. + Args: mask: The current distribution is generated at the edges of the mask, i.e. any points where shifting the mask by one cell in any direction diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py index 7d1f651..1074e2b 100644 --- a/meanas/fdfd/waveguide_2d.py +++ b/meanas/fdfd/waveguide_2d.py @@ -175,8 +175,6 @@ if the result is introduced into a space with a discretized z-axis. """ -# TODO update module docs - from typing import Any from collections.abc import Sequence import numpy @@ -339,7 +337,7 @@ def normalized_fields_e( mu: vfdslice | None = None, prop_phase: float = 0, ) -> tuple[vcfdslice_t, vcfdslice_t]: - """ + r""" Given a vector `e_xy` containing the vectorized E_x and E_y fields, returns normalized, vectorized E and H fields for the system. @@ -357,6 +355,21 @@ def normalized_fields_e( Returns: `(e, h)`, where each field is vectorized, normalized, and contains all three vector components. + + Notes: + `e_xy` is only the transverse electric eigenvector. This helper first + reconstructs the full three-component `E` and `H` fields with `exy2e(...)` + and `exy2h(...)`, then normalizes them to unit forward power using + `_normalized_fields(...)`. + + The normalization target is + + $$ + \Re\left[\mathrm{inner\_product}(e, h, \mathrm{conj\_h}=True)\right] = 1, + $$ + + so the returned fields represent a unit-power forward mode under the + discrete Yee-grid Poynting inner product. """ e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy h = exy2h(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ e_xy @@ -374,7 +387,7 @@ def normalized_fields_h( mu: vfdslice | None = None, prop_phase: float = 0, ) -> tuple[vcfdslice_t, vcfdslice_t]: - """ + r""" Given a vector `h_xy` containing the vectorized H_x and H_y fields, returns normalized, vectorized E and H fields for the system. @@ -392,6 +405,13 @@ def normalized_fields_h( Returns: `(e, h)`, where each field is vectorized, normalized, and contains all three vector components. + + Notes: + This is the `H_x/H_y` analogue of `normalized_fields_e(...)`. The final + normalized mode should describe the same physical solution, but because + the overall complex phase and sign are chosen heuristically, + `normalized_fields_e(...)` and `normalized_fields_h(...)` need not return + identical representatives for nearly symmetric modes. """ e = hxy2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ h_xy h = hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) @ h_xy @@ -409,7 +429,25 @@ def _normalized_fields( mu: vfdslice | None = None, prop_phase: float = 0, ) -> tuple[vcfdslice_t, vcfdslice_t]: - # TODO documentation + r""" + Normalize a reconstructed waveguide mode to unit forward power. + + The eigenproblem solved by `solve_mode(s)` determines only the mode shape and + propagation constant. The overall complex amplitude and sign are still free. + This helper fixes those remaining degrees of freedom in two steps: + + 1. Compute the discrete longitudinal Poynting flux with + `inner_product(e, h, conj_h=True)`, including the half-cell longitudinal + phase adjustment controlled by `prop_phase`. + 2. Multiply both fields by a scalar chosen so that the real forward power is + `1`, then choose a reproducible phase/sign representative by making a + dominant-energy sample real and using a weighted quadrant sum to break + mirror-symmetry ties. + + The sign heuristic is intentionally pragmatic rather than fundamental: it is + only there to make downstream tests and source/overlap construction choose a + consistent representative when the physical mode is symmetric. + """ shape = [s.size for s in dxes[0]] # Find time-averaged Sz and normalize to it @@ -921,7 +959,7 @@ def solve_mode( return vcfdfield2_t(e_xys[0]), wavenumbers[0] -def inner_product( # TODO documentation +def inner_product( e1: vcfdfield2, h2: vcfdfield2, dxes: dx_lists2_t, @@ -929,6 +967,36 @@ def inner_product( # TODO documentation conj_h: bool = False, trapezoid: bool = False, ) -> complex: + r""" + Compute the discrete waveguide overlap / Poynting inner product. + + This is the 2D transverse integral corresponding to the time-averaged + longitudinal Poynting flux, + + $$ + \frac{1}{2}\int (E_x H_y - E_y H_x) \, dx \, dy + $$ + + with the Yee-grid staggering and optional propagation-phase adjustment used + by the waveguide helpers in this module. + + Args: + e1: Vectorized electric field, typically from `exy2e(...)` or + `normalized_fields_e(...)`. + h2: Vectorized magnetic field, typically from `hxy2h(...)`, + `exy2h(...)`, or one of the normalization helpers. + dxes: Two-dimensional Yee-grid spacing lists `[dx_e, dx_h]`. + prop_phase: Phase advance over one propagation cell. This is used to + shift the H field into the same longitudinal reference plane as the + E field. + conj_h: Whether to conjugate `h2` before forming the overlap. Use + `True` for the usual time-averaged power normalization. + trapezoid: Whether to use trapezoidal quadrature instead of the default + rectangular Yee-cell sum. + + Returns: + Complex overlap / longitudinal power integral. + """ shape = [s.size for s in dxes[0]] @@ -951,5 +1019,3 @@ def inner_product( # TODO documentation Sz_b = E1[1] * H2[0] * dxes_real[0][0][:, None] * dxes_real[1][1][None, :] Sz = 0.5 * (Sz_a.sum() - Sz_b.sum()) return Sz - - diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py index 19975db..66cc7cc 100644 --- a/meanas/fdfd/waveguide_3d.py +++ b/meanas/fdfd/waveguide_3d.py @@ -3,6 +3,21 @@ Tools for working with waveguide modes in 3D domains. This module relies heavily on `waveguide_2d` and mostly just transforms its parameters into 2D equivalents and expands the results back into 3D. + +The intended workflow is: + +1. Select a single-cell slice normal to the propagation axis. +2. Solve the corresponding 2D mode problem with `solve_mode(...)`. +3. Turn that mode into a one-sided source with `compute_source(...)`. +4. Build an overlap window with `compute_overlap_e(...)` for port readout. + +`polarity` is part of the public convention throughout this module: + +- `+1` means forward propagation toward increasing index along `axis` +- `-1` means backward propagation toward decreasing index along `axis` + +That same convention controls which side of the selected slice is used for the +overlap window and how the expanded field is phased. """ from typing import Any, cast import warnings @@ -26,9 +41,9 @@ def solve_mode( epsilon: fdfield, mu: fdfield | None = None, ) -> dict[str, complex | NDArray[complexfloating]]: - """ + r""" Given a 3D grid, selects a slice from the grid and attempts to - solve for an eigenmode propagating through that slice. + solve for an eigenmode propagating through that slice. Args: mode_number: Number of the mode, 0-indexed @@ -37,19 +52,22 @@ def solve_mode( axis: Propagation axis (0=x, 1=y, 2=z) polarity: Propagation direction (+1 for +ve, -1 for -ve) slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use - as the waveguide cross-section. `slices[axis]` should select only one item. + as the waveguide cross-section. `slices[axis]` must select exactly one item. epsilon: Dielectric constant mu: Magnetic permeability (default 1 everywhere) Returns: - ``` - { - 'E': NDArray[complexfloating], - 'H': NDArray[complexfloating], - 'wavenumber': complex, - 'wavenumber_2d': complex, - } - ``` + Dictionary containing: + + - `E`: full-grid electric field for the solved mode + - `H`: full-grid magnetic field for the solved mode + - `wavenumber`: propagation constant corrected for the discretized + propagation axis + - `wavenumber_2d`: propagation constant of the reduced 2D eigenproblem + + Notes: + The returned fields are normalized through the `waveguide_2d` + normalization convention before being expanded back to 3D. """ if mu is None: mu = numpy.ones_like(epsilon) @@ -139,7 +157,14 @@ def compute_source( mu: Magnetic permeability (default 1 everywhere) Returns: - J distribution for the unidirectional source + `J` distribution for a one-sided electric-current source. + + Notes: + The source is built from the expanded mode field and a boundary-source + operator. The resulting current is intended to be injected with the + same sign convention used elsewhere in the package: + + `E -= dt * J / epsilon` """ E_expanded = expand_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis, polarity=polarity, slices=slices) @@ -167,10 +192,33 @@ def compute_overlap_e( slices: Sequence[slice], omega: float, ) -> cfdfield_t: - """ - Given an eigenmode obtained by `solve_mode`, calculates an overlap_e for the - mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn) - [assumes reflection symmetry]. + r""" + Build an overlap field for projecting another 3D electric field onto a mode. + + The returned field is intended for the discrete overlap expression + + $$ + \sum \mathrm{overlap\_e} \; E_\mathrm{other}^* + $$ + + where the sum is over the full Yee-grid field storage. + + The construction uses a two-cell window immediately upstream of the selected + slice: + + - for `polarity=+1`, the two cells just before `slices[axis].start` + - for `polarity=-1`, the two cells just after `slices[axis].stop` + + The window is clipped to the simulation domain if necessary. A clipped but + non-empty window raises `RuntimeWarning`; an empty window raises + `ValueError`. + + The derivation below assumes reflection symmetry and the standard waveguide + overlap relation involving + + $$ + \int ((E \times H_\mathrm{mode}) + (E_\mathrm{mode} \times H)) \cdot dn. + $$ E x H_mode + E_mode x H -> Ex Hmy - EyHmx + Emx Hy - Emy Hx (Z-prop) @@ -183,8 +231,6 @@ def compute_overlap_e( B/wu (Ex Emx + Ey Emy) - j/wu (Ex dx Emz + Ey dy Emz) - TODO: add reference - Args: E: E-field of the mode wavenumber: Wavenumber of the mode @@ -195,7 +241,8 @@ def compute_overlap_e( as the waveguide cross-section. slices[axis] should select only one item. Returns: - overlap_e such that `numpy.sum(overlap_e * other_e.conj())` computes the overlap integral + `overlap_e` normalized so that `numpy.sum(overlap_e * E.conj()) == 1` + over the retained overlap window. """ slices = tuple(slices) @@ -240,7 +287,7 @@ def expand_e( polarity: int, slices: Sequence[slice], ) -> cfdfield_t: - """ + r""" Given an eigenmode obtained by `solve_mode`, expands the E-field from the 2D slice where the mode was calculated to the entire domain (along the propagation axis). This assumes the epsilon cross-section remains constant throughout the @@ -258,6 +305,16 @@ def expand_e( Returns: `E`, with the original field expanded along the specified `axis`. + + Notes: + This helper assumes that the waveguide cross-section remains constant + along the propagation axis and applies the phase factor + + $$ + e^{-i \, \mathrm{polarity} \, wavenumber \, \Delta z} + $$ + + to each copied slice. """ slices = tuple(slices) diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index 597e1cb..caedfaf 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -2,63 +2,125 @@ r""" Operators and helper functions for cylindrical waveguides with unchanging cross-section. Waveguide operator is derived according to 10.1364/OL.33.001848. -The curl equations in the complex coordinate system become + +As in `waveguide_2d`, the propagation dependence is separated from the +transverse solve. Here the propagation coordinate is the bend angle `\theta`, +and the fields are assumed to have the form + +$$ +\vec{E}(r, y, \theta), \vec{H}(r, y, \theta) \propto e^{-\imath m \theta}, +$$ + +where `m` is the angular wavenumber returned by `solve_mode(s)`. It is often +convenient to introduce the corresponding linear wavenumber + +$$ +\beta = \frac{m}{r_{\min}}, +$$ + +so that the cylindrical problem resembles the straight-waveguide problem with +additional metric factors. + +Those metric factors live on the staggered radial Yee grids. If the left edge of +the computational window is at `r = r_{\min}`, define the electric-grid and +magnetic-grid radial sample locations by $$ \begin{aligned} --\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta frac{E_y}{\tilde{t}_x} \\ --\imath \omega \mu_{yy} H_y &= -\imath \beta E_x - \frac{1}{\hat{t}_x} \tilde{\partial}_x \tilde{t}_x E_z \\ --\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\ -\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\ -\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\ -\imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ +r_a(n) &= r_{\min} + \sum_{j \le n} \Delta r_{e, j}, \\ +r_b\!\left(n + \tfrac{1}{2}\right) &= r_{\min} + \tfrac{1}{2}\Delta r_{e, n} + + \sum_{j < n} \Delta r_{h, j}, \end{aligned} $$ -where $t_x = 1 + \frac{\Delta_{x, m}}{R_0}$ is the grid spacing adjusted by the nominal radius $R0$. - -Rewrite the last three equations as +and from them the diagonal metric matrices $$ \begin{aligned} -\imath \beta H_y &= \imath \omega \hat{t}_x \epsilon_{xx} E_x - \hat{t}_x \hat{\partial}_y H_z \\ -\imath \beta H_x &= -\imath \omega \hat{t}_x \epsilon_{yy} E_y - \hat{t}_x \hat{\partial}_x H_z \\ -\imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\ +T_a &= \operatorname{diag}(r_a / r_{\min}), \\ +T_b &= \operatorname{diag}(r_b / r_{\min}). \end{aligned} $$ -The derivation then follows the same steps as the straight waveguide, leading to the eigenvalue problem - -$$ -\beta^2 \begin{bmatrix} E_x \\ - E_y \end{bmatrix} = - (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\ - 0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} + - \begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\ - T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \mu_{zz}^{-1} - \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} + - \begin{bmatrix} \tilde{\partial}_x \\ - \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1} - \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix}) - \begin{bmatrix} E_x \\ - E_y \end{bmatrix} -$$ - -which resembles the straight waveguide eigenproblem with additonal $T_a$ and $T_b$ terms. These -are diagonal matrices containing the $t_x$ values: +With the same forward/backward derivative notation used in `waveguide_2d`, the +coordinate-transformed discrete curl equations used here are $$ \begin{aligned} -T_a &= 1 + \frac{\Delta_{x, m }}{R_0} -T_b &= 1 + \frac{\Delta_{x, m + \frac{1}{2} }}{R_0} +-\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\ +-\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r + - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\ +-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \\ +\imath \beta H_y &= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\ +\imath \beta H_r &= \imath \omega T_a \epsilon_{yy} E_y + - T_b T_a^{-1} \hat{\partial}_r (T_b H_z), \\ +\imath \omega E_z &= T_a \epsilon_{zz}^{-1} + \left(\hat{\partial}_r H_y - \hat{\partial}_y H_r\right). \end{aligned} - - -TODO: consider 10.1364/OE.20.021583 for an alternate approach $$ -As in the straight waveguide case, all the functions in this file assume a 2D grid - (i.e. `dxes = [[[dr_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`). +The first three equations are the cylindrical analogue of the straight-guide +relations for `H_r`, `H_y`, and `H_z`. The next two are the metric-weighted +versions of the straight-guide identities for `\imath \beta H_y` and +`\imath \beta H_r`, and the last equation plays the same role as the +longitudinal `E_z` reconstruction in `waveguide_2d`. + +Following the same elimination steps as in `waveguide_2d`, apply +`\imath \beta \tilde{\partial}_r` and `\imath \beta \tilde{\partial}_y` to the +equation for `E_z`, substitute for `\imath \beta H_r` and `\imath \beta H_y`, +and then eliminate `H_z` with + +$$ +H_z = \frac{1}{-\imath \omega \mu_{zz}} +\left(\tilde{\partial}_r E_y - \tilde{\partial}_y E_r\right). +$$ + +This yields the transverse electric eigenproblem implemented by +`cylindrical_operator(...)`: + +$$ +\beta^2 +\begin{bmatrix} E_r \\ E_y \end{bmatrix} += +\left( +\omega^2 +\begin{bmatrix} +T_b^2 \mu_{yy} \epsilon_{xx} & 0 \\ +0 & T_a^2 \mu_{xx} \epsilon_{yy} +\end{bmatrix} ++ +\begin{bmatrix} +-T_b \mu_{yy} \hat{\partial}_y \\ + T_a \mu_{xx} \hat{\partial}_x +\end{bmatrix} +T_b \mu_{zz}^{-1} +\begin{bmatrix} +-\tilde{\partial}_y & \tilde{\partial}_x +\end{bmatrix} ++ +\begin{bmatrix} +\tilde{\partial}_x \\ +\tilde{\partial}_y +\end{bmatrix} +T_a \epsilon_{zz}^{-1} +\begin{bmatrix} +\hat{\partial}_x T_b \epsilon_{xx} & +\hat{\partial}_y T_a \epsilon_{yy} +\end{bmatrix} +\right) +\begin{bmatrix} E_r \\ E_y \end{bmatrix}. +$$ + +Since `\beta = m / r_{\min}`, the solver implemented in this file returns the +angular wavenumber `m`, while the operator itself is most naturally written in +terms of the linear quantity `\beta`. The helpers below reconstruct the full +field components from the solved transverse eigenvector and then normalize the +mode to unit forward power with the same discrete longitudinal Poynting inner +product used by `waveguide_2d`. + +As in the straight-waveguide case, all functions here assume a 2D grid: + +`dxes = [[[dr_e_0, dr_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`. """ from typing import Any, cast from collections.abc import Sequence @@ -94,17 +156,18 @@ def cylindrical_operator( \begin{bmatrix} \tilde{\partial}_x \\ \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1} \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix}) - \begin{bmatrix} E_x \\ + \begin{bmatrix} E_r \\ E_y \end{bmatrix} $$ for use with a field vector of the form `[E_r, E_y]`. This operator can be used to form an eigenvalue problem of the form - A @ [E_r, E_y] = wavenumber**2 * [E_r, E_y] + A @ [E_r, E_y] = beta**2 * [E_r, E_y] which can then be solved for the eigenmodes of the system - (an `exp(-i * wavenumber * theta)` theta-dependence is assumed for the fields). + (an `exp(-i * angular_wavenumber * theta)` theta-dependence is assumed for + the fields, with `beta = angular_wavenumber / rmin`). (NOTE: See module docs and 10.1364/OL.33.001848) @@ -270,7 +333,7 @@ def exy2h( mu: vfdslice | None = None ) -> sparse.sparray: """ - Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, + Operator which transforms the vector `e_xy` containing the vectorized E_r and E_y fields, into a vectorized H containing all three H components Args: @@ -298,11 +361,11 @@ def exy2e( epsilon: vfdslice, ) -> sparse.sparray: """ - Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, + Operator which transforms the vector `e_xy` containing the vectorized E_r and E_y fields, into a vectorized E containing all three E components Unlike the straight waveguide case, the H_z components do not cancel and must be calculated - from E_x and E_y in order to then calculate E_z. + from E_r and E_y in order to then calculate E_z. Args: angular_wavenumber: Wavenumber assuming fields have theta-dependence of @@ -360,9 +423,10 @@ def e2h( This operator is created directly from the initial coordinate-transformed equations: $$ \begin{aligned} - \imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\ - \imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\ - \imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ + -\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\ + -\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r + - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\ + -\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \end{aligned} $$ @@ -397,15 +461,18 @@ def dxes2T( rmin: float, ) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]: r""" - Returns the $T_a$ and $T_b$ diagonal matrices which are used to apply the cylindrical - coordinate transformation in various operators. + Construct the cylindrical metric matrices $T_a$ and $T_b$. + + `T_a` is sampled on the E-grid radial locations, while `T_b` is sampled on + the staggered H-grid radial locations. These are the diagonal matrices that + convert the straight-waveguide algebra into its cylindrical counterpart. Args: dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) rmin: Radius at the left edge of the simulation domain (at minimum 'x') Returns: - Sparse matrix representations of the operators Ta and Tb + Sparse diagonal matrices `(T_a, T_b)`. """ ra = rmin + numpy.cumsum(dxes[0][0]) # Radius at Ey points rb = rmin + dxes[0][0] / 2.0 + numpy.cumsum(dxes[1][0]) # Radius at Ex points @@ -427,12 +494,12 @@ def normalized_fields_e( mu: vfdslice | None = None, prop_phase: float = 0, ) -> tuple[vcfdslice_t, vcfdslice_t]: - """ - Given a vector `e_xy` containing the vectorized E_x and E_y fields, - returns normalized, vectorized E and H fields for the system. + r""" + Given a vector `e_xy` containing the vectorized E_r and E_y fields, + returns normalized, vectorized E and H fields for the system. Args: - e_xy: Vector containing E_x and E_y fields + e_xy: Vector containing E_r and E_y fields angular_wavenumber: Wavenumber assuming fields have theta-dependence of `exp(-i * angular_wavenumber * theta)`. It should satisfy `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy` @@ -447,6 +514,18 @@ def normalized_fields_e( Returns: `(e, h)`, where each field is vectorized, normalized, and contains all three vector components. + + Notes: + The normalization step is delegated to `_normalized_fields(...)`, which + enforces unit forward power under the discrete inner product + + $$ + \frac{1}{2}\int (E_r H_y^* - E_y H_r^*) \, dr \, dy. + $$ + + The angular wavenumber `m` is first converted into the full three-component + fields, then the overall complex phase and sign are fixed so the result is + reproducible for symmetric modes. """ e = exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy h = exy2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, mu=mu) @ e_xy @@ -465,8 +544,30 @@ def _normalized_fields( mu: vfdslice | None = None, prop_phase: float = 0, ) -> tuple[vcfdslice_t, vcfdslice_t]: + r""" + Normalize a cylindrical waveguide mode to unit forward power. + + The cylindrical helpers reuse the straight-waveguide inner product after the + field reconstruction step. The extra metric factors have already been folded + into the reconstructed `e`/`h` fields through `dxes2T(...)` and the + cylindrical `exy2e(...)` / `exy2h(...)` operators, so the same discrete + longitudinal Poynting integral can be used here. + + The normalization procedure is: + + 1. Flip the magnetic field sign so the reconstructed `(e, h)` pair follows + the same forward-power convention as `waveguide_2d`. + 2. Compute the time-averaged forward power with + `waveguide_2d.inner_product(..., conj_h=True)`. + 3. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power. + 4. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic + so symmetric modes choose a stable representative. + + `prop_phase` has the same meaning as in `waveguide_2d`: it compensates for + the half-cell longitudinal staggering between the E and H fields when the + propagation direction is itself discretized. + """ h *= -1 - # TODO documentation for normalized_fields shape = [s.size for s in dxes[0]] # Find time-averaged Sz and normalize to it diff --git a/meanas/fdtd/__init__.py b/meanas/fdtd/__init__.py index 2334338..2b15c59 100644 --- a/meanas/fdtd/__init__.py +++ b/meanas/fdtd/__init__.py @@ -168,7 +168,44 @@ t=0 (assuming the source is off for t<0 this gives $\sim 10^{-3}$ error at t=0). Boundary conditions =================== -# TODO notes about boundaries / PMLs + +`meanas.fdtd` exposes two boundary-related building blocks: + +- `conducting_boundary(...)` for simple perfect-electric-conductor style field + clamping at one face of the domain. +- `cpml_params(...)` and `updates_with_cpml(...)` for convolutional perfectly + matched layers (CPMLs) attached to one or more faces of the Yee grid. + +`updates_with_cpml(...)` accepts a three-by-two table of CPML parameter blocks: + +``` +cpml_params[axis][polarity_index] +``` + +where `axis` is `0`, `1`, or `2` and `polarity_index` corresponds to `(-1, +1)`. +Passing `None` for one entry disables CPML on that face while leaving the other +directions unchanged. This is how mixed boundary setups such as "absorbing in x, +periodic in y/z" are expressed. + +When comparing an FDTD run against an FDFD solve, use the same stretched +coordinate system in both places: + +1. Build the FDTD update with the desired CPML parameters. +2. Stretch the FDFD `dxes` with the matching SCPML transform. +3. Compare the extracted phasor against the FDFD field or residual on those + stretched `dxes`. + +The electric-current sign convention used throughout the examples and tests is + +$$ +E \leftarrow E - \Delta_t J / \epsilon +$$ + +which matches the FDFD right-hand side + +$$ +-i \omega J. +$$ """ from .base import ( diff --git a/meanas/fdtd/energy.py b/meanas/fdtd/energy.py index 76888ca..6df30dc 100644 --- a/meanas/fdtd/energy.py +++ b/meanas/fdtd/energy.py @@ -4,7 +4,21 @@ from ..fdmath import dx_lists_t, fdfield_t, fdfield from ..fdmath.functional import deriv_back -# TODO documentation +""" +Energy- and flux-accounting helpers for Yee-grid FDTD fields. + +These functions complement the derivation in `meanas.fdtd`: + +- `poynting(...)` and `poynting_divergence(...)` evaluate the discrete flux terms + from the exact time-domain Poynting identity. +- `energy_hstep(...)` / `energy_estep(...)` evaluate the two staggered energy + expressions. +- `delta_energy_*` helpers evaluate the corresponding energy changes between + adjacent half-steps. + +The helpers are intended for diagnostics, regression tests, and consistency +checks between source work, field energy, and flux through cell faces. +""" def poynting( @@ -252,13 +266,23 @@ def delta_energy_j( e1: fdfield, dxes: dx_lists_t | None = None, ) -> fdfield_t: - """ - Calculate + r""" + Calculate the electric-current work term $J \cdot E$ on the Yee grid. - 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$). + This is the source contribution that appears beside the flux divergence in + the discrete Poynting identities documented in `meanas.fdtd`. + Note that each value of `J` contributes twice in a full Yee cycle (once per + half-step energy balance) even though it directly changes `E` only once. + Args: + j0: Electric-current density sampled at the same half-step as the + current work term. + e1: Electric field sampled at the matching integer timestep. + dxes: Grid description; defaults to unit spacing. + + Returns: + Per-cell source-work contribution with the scalar field shape. """ if dxes is None: dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2)) @@ -277,6 +301,20 @@ def dxmul( mu: fdfield | float | None = None, dxes: dx_lists_t | None = None, ) -> fdfield_t: + """ + Multiply E- and H-like field products by material weights and cell volumes. + + Args: + ee: Three-component electric-field product, such as `e0 * e2`. + hh: Three-component magnetic-field product, such as `h1 * h1`. + epsilon: Electric material weight; defaults to `1`. + mu: Magnetic material weight; defaults to `1`. + dxes: Grid description; defaults to unit spacing. + + Returns: + Scalar field containing the weighted electric plus magnetic contribution + for each Yee cell. + """ if epsilon is None: epsilon = 1 if mu is None: diff --git a/meanas/fdtd/pml.py b/meanas/fdtd/pml.py index dec3d83..bf61b4e 100644 --- a/meanas/fdtd/pml.py +++ b/meanas/fdtd/pml.py @@ -1,9 +1,19 @@ """ -PML implementations +Convolutional perfectly matched layer (CPML) support for FDTD updates. -#TODO discussion of PMLs -#TODO cpml documentation +The helpers in this module construct per-face CPML parameters and then wrap the +standard Yee updates with the additional auxiliary `psi` fields needed by the +CPML recurrence. +The intended call pattern is: + +1. Build a `cpml_params[axis][polarity_index]` table with `cpml_params(...)`. +2. Pass that table into `updates_with_cpml(...)` together with `dt`, `dxes`, and + `epsilon`. +3. Advance the returned `update_E` / `update_H` closures in the simulation loop. + +Each face can be enabled or disabled independently by replacing one table entry +with `None`. """ # TODO retest pmls! @@ -32,6 +42,29 @@ def cpml_params( ma: float = 1, cfs_alpha: float = 0, ) -> dict[str, Any]: + """ + Construct the parameter block for one CPML face. + + Args: + axis: Which Cartesian axis the CPML is normal to (`0`, `1`, or `2`). + polarity: Which face along that axis (`-1` for the low-index face, + `+1` for the high-index face). + dt: Timestep used by the Yee update. + thickness: Number of Yee cells occupied by the CPML region. + ln_R_per_layer: Logarithmic attenuation target per layer. + epsilon_eff: Effective permittivity used when choosing the CPML scaling. + mu_eff: Effective permeability used when choosing the CPML scaling. + m: Polynomial grading exponent for `sigma` and `kappa`. + ma: Polynomial grading exponent for the complex-frequency shift `alpha`. + cfs_alpha: Maximum complex-frequency shift parameter. + + Returns: + Dictionary with: + + - `param_e`: `(p0, p1, p2)` arrays for the E update + - `param_h`: `(p0, p1, p2)` arrays for the H update + - `region`: slice tuple selecting the CPML cells on that face + """ if axis not in range(3): raise Exception(f'Invalid axis: {axis}') @@ -98,6 +131,27 @@ def updates_with_cpml( dtype: DTypeLike = numpy.float32, ) -> tuple[Callable[[fdfield_t, fdfield_t, fdfield_t], None], Callable[[fdfield_t, fdfield_t, fdfield_t], None]]: + """ + Build Yee-step update closures augmented with CPML terms. + + Args: + cpml_params: Three-by-two sequence indexed as `[axis][polarity_index]`. + Entries are the dictionaries returned by `cpml_params(...)`; use + `None` to disable CPML on one face. + dt: Timestep. + dxes: Yee-grid spacing lists `[dx_e, dx_h]`. + epsilon: Electric material distribution used by the E update. + dtype: Storage dtype for the auxiliary CPML state arrays. + + Returns: + `(update_E, update_H)` closures with the same call shape as the basic + Yee updates: + + - `update_E(e, h, epsilon)` + - `update_H(e, h, mu)` + + The closures retain the CPML auxiliary state internally. + """ Dfx, Dfy, Dfz = deriv_forward(dxes[1]) Dbx, Dby, Dbz = deriv_back(dxes[1])