[eme / examples] add EME examples

README.md

@@ -163,6 +163,11 @@ The tracked examples under `examples/` are the intended entry points for users:
   guide, with late-time monitor slices, guided-core windows, and mode-weighted
   errors compared directly against real fields reconstructed from the matching
   FDFD solution, plus a guided-mode / orthogonal-residual split.
+- `examples/eme.py`: straight-interface mode matching / EME, including port
+  mode solving, interface scattering, and modal field visualization.
+- `examples/eme_bend.py`: straight-to-bent waveguide mode matching with
+  cylindrical bend modes, interface scattering, and a cascaded bend-network
+  example built with `scikit-rf`.
 - `examples/fdfd.py`: direct frequency-domain waveguide excitation and overlap /
   Poynting analysis without a time-domain run.
 

docs/index.md

@@ -24,6 +24,10 @@ Relevant starting examples:
 - `examples/waveguide_real.py` for real-valued continuous-wave FDTD compared
   against real fields reconstructed from an FDFD solution, including guided-core,
   mode-weighted, and guided-mode / residual comparisons
+- `examples/eme.py` for straight-interface mode matching / EME and modal
+  scattering between two nearby waveguide cross-sections
+- `examples/eme_bend.py` for straight-to-bent mode matching with cylindrical
+  bend modes and a cascaded bend-network example
 - `examples/fdfd.py` for direct frequency-domain waveguide excitation
 
 For solver equivalence, prefer the phasor-based examples first. They compare

docs/stylesheets/extra.css

@@ -11,3 +11,42 @@
 .md-typeset h3 code {
   word-break: break-word;
 }
+
+[data-md-color-scheme="slate"] {
+  --md-default-bg-color: #0f141c;
+  --md-default-fg-color: #e8eef7;
+  --md-default-fg-color--light: #b3bfd1;
+  --md-default-fg-color--lighter: #7f8ba0;
+  --md-default-fg-color--lightest: #5d6880;
+  --md-code-bg-color: #111923;
+  --md-code-fg-color: #e4edf8;
+  --md-accent-fg-color: #7dd3fc;
+}
+
+[data-md-color-scheme="slate"] .md-header,
+[data-md-color-scheme="slate"] .md-tabs {
+  background: linear-gradient(90deg, #111923 0%, #162235 100%);
+}
+
+[data-md-color-scheme="slate"] .md-typeset pre > code,
+[data-md-color-scheme="slate"] .md-typeset code {
+  border: 1px solid rgba(125, 211, 252, 0.14);
+}
+
+[data-md-color-scheme="slate"] .md-typeset table:not([class]) {
+  background: rgba(255, 255, 255, 0.015);
+}
+
+[data-md-color-scheme="slate"] .md-typeset table:not([class]) th {
+  background: rgba(125, 211, 252, 0.08);
+}
+
+[data-md-color-scheme="slate"] .md-typeset .admonition,
+[data-md-color-scheme="slate"] .md-typeset details {
+  background: rgba(255, 255, 255, 0.02);
+  border-color: rgba(125, 211, 252, 0.2);
+}
+
+[data-md-color-scheme="slate"] .md-typeset .arithmatex {
+  padding: 0.1rem 0;
+}

examples/eme.py

@@ -0,0 +1,217 @@
+"""
+Mode-matching / EME example for a straight rib-waveguide interface.
+
+This example shows the intended user-facing workflow for `meanas.fdfd.eme` on a
+simple straight interface:
+
+1. build two nearby waveguide cross-sections on a Yee grid,
+2. solve a small set of guided modes on each side,
+3. normalize those modes into E/H port fields,
+4. assemble the interface scattering matrix with `meanas.fdfd.eme.get_s(...)`,
+5. inspect the dominant modal coupling numerically and visually.
+"""
+
+from __future__ import annotations
+
+import importlib
+
+import numpy
+from numpy import pi
+
+import gridlock
+from gridlock import Extent
+
+from meanas.fdfd import eme, waveguide_2d
+from meanas.fdmath import unvec
+
+
+WL = 1310.0
+DX = 40.0
+WIDTH = 400.0
+THF = 161.0
+THP = 77.0
+EPS_SI = 3.51 ** 2
+EPS_OX = 1.453 ** 2
+MODE_NUMBERS = numpy.array([0])
+
+
+def require_optional(name: str, package_name: str | None = None):
+    package_name = package_name or name
+    try:
+        return importlib.import_module(name)
+    except ImportError as exc:  # pragma: no cover - user environment guard
+        raise SystemExit(
+            f"This example requires the optional package '{package_name}'. "
+            "Install example dependencies with `pip install -e './meanas[examples]'`.",
+        ) from exc
+
+
+def build_geometry(
+        *,
+        dx: float = DX,
+        width: float = WIDTH,
+        thf: float = THF,
+        thp: float = THP,
+        eps_si: float = EPS_SI,
+        eps_ox: float = EPS_OX,
+        ) -> tuple[gridlock.Grid, numpy.ndarray, list[list[numpy.ndarray]], float]:
+    x0 = (width / 2) % dx
+    omega = 2 * pi / WL
+
+    grid = gridlock.Grid(
+        [
+            numpy.arange(-800, 800 + dx, dx),
+            numpy.arange(-400, 400 + dx, dx),
+            numpy.arange(-2 * dx, 2 * dx + dx, dx),
+        ],
+        periodic=True,
+    )
+    epsilon = grid.allocate(eps_ox)
+
+    grid.draw_cuboid(
+        epsilon,
+        foreground=eps_si,
+        x=Extent(center=x0, span=width + 1200),
+        y=Extent(min=0, max=thf),
+        z=Extent(min=-1e6, max=0),
+    )
+    grid.draw_cuboid(
+        epsilon,
+        foreground=eps_ox,
+        x=Extent(max=-width / 2, span=300),
+        y=Extent(min=thp, max=1e6),
+        z=Extent(min=-1e6, max=0),
+    )
+    grid.draw_cuboid(
+        epsilon,
+        foreground=eps_ox,
+        x=Extent(min=width / 2, span=300),
+        y=Extent(min=thp, max=1e6),
+        z=Extent(min=-1e6, max=0),
+    )
+
+    grid.draw_cuboid(
+        epsilon,
+        foreground=eps_si,
+        x=Extent(max=-(width / 2 + 600), span=240),
+        y=Extent(min=0, max=thf),
+        z=Extent(min=0, max=1e6),
+    )
+    grid.draw_cuboid(
+        epsilon,
+        foreground=eps_si,
+        x=Extent(max=width / 2 + 600, span=240),
+        y=Extent(min=0, max=thf),
+        z=Extent(min=0, max=1e6),
+    )
+
+    dxes = [grid.dxyz, grid.autoshifted_dxyz()]
+    dxes_2d = [[d[0], d[1]] for d in dxes]
+    return grid, epsilon, dxes_2d, omega
+
+
+def solve_cross_section_modes(
+        epsilon_slice: numpy.ndarray,
+        *,
+        omega: float,
+        dxes_2d: list[list[numpy.ndarray]],
+        mode_numbers: numpy.ndarray = MODE_NUMBERS,
+        ) -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], numpy.ndarray]:
+    e_xys, wavenumbers = waveguide_2d.solve_modes(
+        epsilon=epsilon_slice.ravel(),
+        omega=omega,
+        dxes=dxes_2d,
+        mode_numbers=mode_numbers,
+    )
+    eh_fields = [
+        waveguide_2d.normalized_fields_e(
+            e_xy,
+            wavenumber=wavenumber,
+            dxes=dxes_2d,
+            omega=omega,
+            epsilon=epsilon_slice.ravel(),
+        )
+        for e_xy, wavenumber in zip(e_xys, wavenumbers, strict=True)
+    ]
+    return eh_fields, wavenumbers
+
+
+def print_summary(ss: numpy.ndarray, wavenumbers_left: numpy.ndarray, wavenumbers_right: numpy.ndarray, omega: float) -> None:
+    n_left = len(wavenumbers_left)
+    left_neff = numpy.real(wavenumbers_left / omega)
+    right_neff = numpy.real(wavenumbers_right / omega)
+
+    print('left effective indices:', ', '.join(f'{value:.5f}' for value in left_neff[:4]))
+    print('right effective indices:', ', '.join(f'{value:.5f}' for value in right_neff[:4]))
+
+    reflection = abs(ss[0, 0]) ** 2
+    transmission = abs(ss[n_left, 0]) ** 2
+    total_output = numpy.sum(abs(ss[:, 0]) ** 2)
+    print(f'fundamental left-incident reflection |S_00|^2      = {reflection:.6f}')
+    print(f'fundamental left-to-right transmission |S_{n_left},0|^2 = {transmission:.6f}')
+    print(f'fundamental left-incident total output power      = {total_output:.6f}')
+
+    strongest = numpy.argsort(abs(ss[n_left:, 0]) ** 2)[::-1][:3]
+    print('dominant transmitted right-side modes for left mode 0:')
+    for mode_index in strongest:
+        print(f'  mode {mode_index}: |S|^2 = {abs(ss[n_left + mode_index, 0]) ** 2:.6f}')
+
+
+def plot_results(
+        *,
+        pyplot,
+        ss: numpy.ndarray,
+        left_mode: tuple[numpy.ndarray, numpy.ndarray],
+        right_mode: tuple[numpy.ndarray, numpy.ndarray],
+        shape_2d: tuple[int, int],
+        ) -> None:
+    fig_s, ax_s = pyplot.subplots()
+    image = ax_s.imshow(abs(ss) ** 2, origin='lower', cmap='magma')
+    fig_s.colorbar(image, ax=ax_s)
+    ax_s.set_title('Interface scattering magnitude |S|^2')
+    ax_s.set_xlabel('Incident mode index')
+    ax_s.set_ylabel('Outgoing mode index')
+
+    e_left = unvec(left_mode[0], shape_2d)
+    e_right = unvec(right_mode[0], shape_2d)
+    left_intensity = numpy.sum(abs(e_left) ** 2, axis=0).T
+    right_intensity = numpy.sum(abs(e_right) ** 2, axis=0).T
+
+    fig_modes, axes = pyplot.subplots(1, 2, figsize=(10, 4))
+    left_plot = axes[0].imshow(left_intensity, origin='lower', cmap='viridis')
+    fig_modes.colorbar(left_plot, ax=axes[0])
+    axes[0].set_title('Left fundamental mode |E|^2')
+    right_plot = axes[1].imshow(right_intensity, origin='lower', cmap='viridis')
+    fig_modes.colorbar(right_plot, ax=axes[1])
+    axes[1].set_title('Right fundamental mode |E|^2')
+    if pyplot.get_backend().lower().endswith('agg'):
+        pyplot.close(fig_s)
+        pyplot.close(fig_modes)
+    else:
+        pyplot.show()
+
+
+def main() -> None:
+    pyplot = require_optional('matplotlib.pyplot', package_name='matplotlib')
+
+    grid, epsilon, dxes_2d, omega = build_geometry()
+    left_slice = epsilon[:, :, :, 1]
+    right_slice = epsilon[:, :, :, -2]
+
+    left_modes, wavenumbers_left = solve_cross_section_modes(left_slice, omega=omega, dxes_2d=dxes_2d)
+    right_modes, wavenumbers_right = solve_cross_section_modes(right_slice, omega=omega, dxes_2d=dxes_2d)
+
+    ss = eme.get_s(left_modes, wavenumbers_left, right_modes, wavenumbers_right, dxes=dxes_2d)
+
+    print_summary(ss, wavenumbers_left, wavenumbers_right, omega)
+    plot_results(
+        pyplot=pyplot,
+        ss=ss,
+        left_mode=left_modes[0],
+        right_mode=right_modes[0],
+        shape_2d=grid.shape[:2],
+    )
+
+
+if __name__ == '__main__':
+    main()

examples/eme_bend.py

@@ -0,0 +1,310 @@
+"""
+Mode-matching / EME example for a straight-to-bent waveguide interface.
+
+This example demonstrates a cylindrical-waveguide EME workflow:
+
+1. build a rib-waveguide cross-section,
+2. solve straight port modes with `waveguide_2d`,
+3. solve bend modes with `waveguide_cyl`,
+4. assemble the straight-to-bend interface scattering matrix with
+   `meanas.fdfd.eme.get_s(...)`,
+5. optionally cascade a straight section, bend section, and interface pair into
+   a compact multiport network using `scikit-rf`.
+"""
+
+from __future__ import annotations
+
+import importlib
+
+import numpy
+from numpy import pi
+from scipy import sparse
+
+import gridlock
+from gridlock import Extent
+
+from meanas.fdfd import eme, waveguide_2d, waveguide_cyl
+from meanas.fdmath import unvec
+
+
+WL = 1310.0
+DX = 40.0
+RADIUS = 6e3
+WIDTH = 400.0
+THF = 161.0
+THP = 77.0
+EPS_SI = 3.51 ** 2
+EPS_OX = 1.453 ** 2
+MODE_NUMBERS = numpy.array([0])
+STRAIGHT_SECTION_LENGTH = 12e3
+BEND_ANGLE = pi / 2
+
+
+def require_optional(name: str, package_name: str | None = None):
+    package_name = package_name or name
+    try:
+        return importlib.import_module(name)
+    except ImportError as exc:  # pragma: no cover - user environment guard
+        raise SystemExit(
+            f"This example requires the optional package '{package_name}'. "
+            "Install example dependencies with `pip install -e './meanas[examples]'`.",
+        ) from exc
+
+
+def build_geometry(
+        *,
+        dx: float = DX,
+        width: float = WIDTH,
+        thf: float = THF,
+        thp: float = THP,
+        eps_si: float = EPS_SI,
+        eps_ox: float = EPS_OX,
+        ) -> tuple[gridlock.Grid, numpy.ndarray, list[list[numpy.ndarray]], float]:
+    x0 = (width / 2) % dx
+    omega = 2 * pi / WL
+
+    grid = gridlock.Grid(
+        [
+            numpy.arange(-800, 800 + dx, dx),
+            numpy.arange(-400, 400 + dx, dx),
+            numpy.arange(-2 * dx, 2 * dx + dx, dx),
+        ],
+        periodic=True,
+    )
+    epsilon = grid.allocate(eps_ox)
+
+    grid.draw_cuboid(
+        epsilon,
+        foreground=eps_si,
+        x=Extent(center=x0, span=width + 1200),
+        y=Extent(min=0, max=thf),
+        z=Extent(min=-1e6, max=0),
+    )
+    grid.draw_cuboid(
+        epsilon,
+        foreground=eps_ox,
+        x=Extent(max=-width / 2, span=300),
+        y=Extent(min=thp, max=1e6),
+        z=Extent(min=-1e6, center=0),
+    )
+    grid.draw_cuboid(
+        epsilon,
+        foreground=eps_ox,
+        x=Extent(min=width / 2, span=300),
+        y=Extent(min=thp, max=1e6),
+        z=Extent(min=-1e6, center=0),
+    )
+
+    dxes = [grid.dxyz, grid.autoshifted_dxyz()]
+    dxes_2d = [[d[0], d[1]] for d in dxes]
+    return grid, epsilon, dxes_2d, omega
+
+
+def solve_straight_modes(
+        epsilon_slice: numpy.ndarray,
+        *,
+        omega: float,
+        dxes_2d: list[list[numpy.ndarray]],
+        mode_numbers: numpy.ndarray = MODE_NUMBERS,
+        ) -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], numpy.ndarray]:
+    e_xys, wavenumbers = waveguide_2d.solve_modes(
+        epsilon=epsilon_slice.ravel(),
+        omega=omega,
+        dxes=dxes_2d,
+        mode_numbers=mode_numbers,
+    )
+    eh_fields = [
+        waveguide_2d.normalized_fields_e(
+            e_xy,
+            wavenumber=wavenumber,
+            dxes=dxes_2d,
+            omega=omega,
+            epsilon=epsilon_slice.ravel(),
+        )
+        for e_xy, wavenumber in zip(e_xys, wavenumbers, strict=True)
+    ]
+    return eh_fields, wavenumbers
+
+
+def solve_bend_modes(
+        epsilon_slice: numpy.ndarray,
+        *,
+        omega: float,
+        dxes_2d: list[list[numpy.ndarray]],
+        rmin: float,
+        mode_numbers: numpy.ndarray = MODE_NUMBERS,
+        ) -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], numpy.ndarray, numpy.ndarray]:
+    e_xys, angular_wavenumbers = waveguide_cyl.solve_modes(
+        epsilon=epsilon_slice.ravel(),
+        omega=omega,
+        dxes=dxes_2d,
+        mode_numbers=mode_numbers,
+        rmin=rmin,
+    )
+    linear_wavenumbers = waveguide_cyl.linear_wavenumbers(
+        e_xys=e_xys,
+        angular_wavenumbers=angular_wavenumbers,
+        rmin=rmin,
+        epsilon=epsilon_slice.ravel(),
+        dxes=dxes_2d,
+    )
+    eh_fields = [
+        waveguide_cyl.normalized_fields_e(
+            e_xy,
+            angular_wavenumber=angular_wavenumber,
+            dxes=dxes_2d,
+            omega=omega,
+            epsilon=epsilon_slice.ravel(),
+            rmin=rmin,
+        )
+        for e_xy, angular_wavenumber in zip(e_xys, angular_wavenumbers, strict=True)
+    ]
+    return eh_fields, linear_wavenumbers, angular_wavenumbers
+
+
+def build_cascaded_network(
+        skrf,
+        *,
+        interface_s: numpy.ndarray,
+        straight_wavenumbers: numpy.ndarray,
+        bend_angular_wavenumbers: numpy.ndarray,
+        n_modes: int,
+        ) -> object:
+    net_sb = skrf.Network(f=[1 / WL], s=interface_s[numpy.newaxis, ...])
+    net_bs = net_sb.copy()
+    net_bs.renumber(numpy.arange(2 * n_modes), numpy.roll(numpy.arange(2 * n_modes), n_modes))
+
+    straight_phase = sparse.diags_array(numpy.exp(-1j * straight_wavenumbers[:n_modes] * STRAIGHT_SECTION_LENGTH))
+    bend_phase = sparse.diags_array(numpy.exp(-1j * bend_angular_wavenumbers[:n_modes] * BEND_ANGLE))
+    zero = numpy.zeros((n_modes, n_modes), dtype=complex)
+
+    straight_s = numpy.block([[zero, straight_phase.toarray()], [straight_phase.toarray(), zero]])
+    bend_s = numpy.block([[zero, bend_phase.toarray()], [bend_phase.toarray(), zero]])
+    net_straight = skrf.Network(f=[1 / WL], s=straight_s[numpy.newaxis, ...])
+    net_bend = skrf.Network(f=[1 / WL], s=bend_s[numpy.newaxis, ...])
+
+    return skrf.network.cascade_list([net_straight, net_sb, net_bend, net_bs, net_straight])
+
+
+def print_summary(
+        interface_s: numpy.ndarray,
+        cascaded_s: numpy.ndarray,
+        straight_wavenumbers: numpy.ndarray,
+        bend_linear_wavenumbers: numpy.ndarray,
+        bend_angular_wavenumbers: numpy.ndarray,
+        omega: float,
+        n_modes: int,
+        ) -> None:
+    straight_neff = numpy.real(straight_wavenumbers / omega)
+    bend_neff = numpy.real(bend_linear_wavenumbers / omega)
+    print('straight effective indices:', ', '.join(f'{value:.5f}' for value in straight_neff[:4]))
+    print('bend effective indices   :', ', '.join(f'{value:.5f}' for value in bend_neff[:4]))
+    print('bend angular wavenumbers :', ', '.join(f'{value:.5e}' for value in bend_angular_wavenumbers[:4]))
+
+    interface_transmission = abs(interface_s[n_modes, 0]) ** 2
+    interface_reflection = abs(interface_s[0, 0]) ** 2
+    print(f'interface fundamental transmission |S_{n_modes},0|^2 = {interface_transmission:.6f}')
+    print(f'interface fundamental reflection   |S_00|^2          = {interface_reflection:.6f}')
+
+    total_cascaded_output = numpy.sum(abs(cascaded_s[:, 0]) ** 2)
+    bend_through = abs(cascaded_s[n_modes, 0]) ** 2
+    bend_reflection = abs(cascaded_s[0, 0]) ** 2
+    print(f'cascaded bend through power        |S_{n_modes},0|^2 = {bend_through:.6f}')
+    print(f'cascaded bend reflection           |S_00|^2          = {bend_reflection:.6f}')
+    print(f'cascaded left-incident total output power            = {total_cascaded_output:.6f}')
+
+
+def plot_results(
+        *,
+        pyplot,
+        interface_s: numpy.ndarray,
+        cascaded_s: numpy.ndarray,
+        straight_mode: tuple[numpy.ndarray, numpy.ndarray],
+        bend_mode: tuple[numpy.ndarray, numpy.ndarray],
+        shape_2d: tuple[int, int],
+        ) -> None:
+    fig_s, axes = pyplot.subplots(1, 2, figsize=(12, 4))
+    interface_plot = axes[0].imshow(abs(interface_s) ** 2, origin='lower', cmap='magma')
+    fig_s.colorbar(interface_plot, ax=axes[0])
+    axes[0].set_title('Straight-to-bend interface |S|^2')
+    axes[0].set_xlabel('Incident mode index')
+    axes[0].set_ylabel('Outgoing mode index')
+
+    cascaded_plot = axes[1].imshow(abs(cascaded_s) ** 2, origin='lower', cmap='magma')
+    fig_s.colorbar(cascaded_plot, ax=axes[1])
+    axes[1].set_title('Cascaded bend network |S|^2')
+    axes[1].set_xlabel('Incident mode index')
+    axes[1].set_ylabel('Outgoing mode index')
+
+    straight_e = unvec(straight_mode[0], shape_2d)
+    bend_e = unvec(bend_mode[0], shape_2d)
+    straight_intensity = numpy.sum(abs(straight_e) ** 2, axis=0).T
+    bend_intensity = numpy.sum(abs(bend_e) ** 2, axis=0).T
+
+    fig_modes, axes_modes = pyplot.subplots(1, 2, figsize=(10, 4))
+    straight_plot = axes_modes[0].imshow(straight_intensity, origin='lower', cmap='viridis')
+    fig_modes.colorbar(straight_plot, ax=axes_modes[0])
+    axes_modes[0].set_title('Straight fundamental mode |E|^2')
+    bend_plot = axes_modes[1].imshow(bend_intensity, origin='lower', cmap='viridis')
+    fig_modes.colorbar(bend_plot, ax=axes_modes[1])
+    axes_modes[1].set_title('Bent fundamental mode |E|^2')
+    if pyplot.get_backend().lower().endswith('agg'):
+        pyplot.close(fig_s)
+        pyplot.close(fig_modes)
+    else:
+        pyplot.show()
+
+
+def main() -> None:
+    pyplot = require_optional('matplotlib.pyplot', package_name='matplotlib')
+    skrf = require_optional('skrf', package_name='scikit-rf')
+
+    grid, epsilon, dxes_2d, omega = build_geometry()
+    epsilon_slice = epsilon[:, :, :, 2]
+    rmin = RADIUS + grid.xyz[0].min()
+
+    straight_modes, straight_wavenumbers = solve_straight_modes(epsilon_slice, omega=omega, dxes_2d=dxes_2d)
+    bend_modes, bend_linear_wavenumbers, bend_angular_wavenumbers = solve_bend_modes(
+        epsilon_slice,
+        omega=omega,
+        dxes_2d=dxes_2d,
+        rmin=rmin,
+    )
+
+    interface_s = eme.get_s(
+        straight_modes,
+        straight_wavenumbers,
+        bend_modes,
+        bend_linear_wavenumbers,
+        dxes=dxes_2d,
+    )
+    cascaded = build_cascaded_network(
+        skrf,
+        interface_s=interface_s,
+        straight_wavenumbers=straight_wavenumbers,
+        bend_angular_wavenumbers=bend_angular_wavenumbers,
+        n_modes=len(MODE_NUMBERS),
+    )
+    cascaded_s = cascaded.s[0]
+
+    print_summary(
+        interface_s,
+        cascaded_s,
+        straight_wavenumbers,
+        bend_linear_wavenumbers,
+        bend_angular_wavenumbers,
+        omega,
+        len(MODE_NUMBERS),
+    )
+    plot_results(
+        pyplot=pyplot,
+        interface_s=interface_s,
+        cascaded_s=cascaded_s,
+        straight_mode=straight_modes[0],
+        bend_mode=bend_modes[0],
+        shape_2d=grid.shape[:2],
+    )
+
+
+if __name__ == '__main__':
+    main()

meanas/fdfd/operators.py

@@ -64,11 +64,11 @@ def e_full(
         epsilon: Vectorized dielectric constant
         mu: Vectorized magnetic permeability (default 1 everywhere).
         pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted
-          as containing a perfect electrical conductor (PEC).
+            as containing a perfect electrical conductor (PEC).
-          The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
+            The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
         pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted
-          as containing a perfect magnetic conductor (PMC).
+            as containing a perfect magnetic conductor (PMC).
-          The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
+            The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
 
     Returns:
         Sparse matrix containing the wave operator.
@@ -148,11 +148,11 @@ def h_full(
         epsilon: Vectorized dielectric constant
         mu: Vectorized magnetic permeability (default 1 everywhere)
         pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted
-           as containing a perfect electrical conductor (PEC).
+            as containing a perfect electrical conductor (PEC).
-           The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
+            The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
         pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted
-           as containing a perfect magnetic conductor (PMC).
+            as containing a perfect magnetic conductor (PMC).
-           The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
+            The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
 
     Returns:
         Sparse matrix containing the wave operator.
@@ -217,11 +217,11 @@ def eh_full(
         epsilon: Vectorized dielectric constant
         mu: Vectorized magnetic permeability (default 1 everywhere)
         pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted
-          as containing a perfect electrical conductor (PEC).
+            as containing a perfect electrical conductor (PEC).
-          The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
+            The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
         pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted
-          as containing a perfect magnetic conductor (PMC).
+            as containing a perfect magnetic conductor (PMC).
-          The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
+            The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
 
     Returns:
         Sparse matrix containing the wave operator.
@@ -267,8 +267,8 @@ def e2h(
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types`
         mu: Vectorized magnetic permeability (default 1 everywhere)
         pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted
-          as containing a perfect magnetic conductor (PMC).
+            as containing a perfect magnetic conductor (PMC).
-          The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
+            The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
 
     Returns:
         Sparse matrix for converting E to H.
@@ -483,4 +483,3 @@ def e_boundary_source(
 #             (numpy.roll(mask, +1, axis=2) != mask))
 
     return sparse.diags_array(jmask.astype(int)) @ full
-

meanas/fdfd/waveguide_3d.py

@@ -52,7 +52,7 @@ 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]` must select exactly one item.
+            as the waveguide cross-section. `slices[axis]` must select exactly one item.
         epsilon: Dielectric constant
         mu: Magnetic permeability (default 1 everywhere)
 
@@ -62,7 +62,7 @@ def solve_mode(
         - `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
+            propagation axis
         - `wavenumber_2d`: propagation constant of the reduced 2D eigenproblem
 
     Notes:

meanas/fdfd/waveguide_cyl.py

@@ -216,13 +216,13 @@ def solve_modes(
      of the bent waveguide with the specified mode number.
 
     Args:
-        mode_number: Number of the mode, 0-indexed
+        mode_numbers: Mode numbers to solve, 0-indexed.
         omega: Angular frequency of the simulation
         dxes: Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types.
-              The first coordinate is assumed to be r, the second is y.
+            The first coordinate is assumed to be r, the second is y.
         epsilon: Dielectric constant
         rmin: Radius of curvature for the simulation. This should be the minimum value of
-               r within the simulation domain.
+            r within the simulation domain.
 
     Returns:
         e_xys: NDArray of vfdfield_t specifying fields. First dimension is mode number.

meanas/fdmath/operators.py

@@ -158,7 +158,7 @@ def cross(
 
     Args:
         B: List `[Bx, By, Bz]` of sparse matrices corresponding to the x, y, z
-           portions of the operator on the left side of the cross product.
+            portions of the operator on the left side of the cross product.
 
     Returns:
         Sparse matrix corresponding to (B x), where x is the cross product.

meanas/fdmath/vectorization.py

@@ -58,7 +58,7 @@ def vec(
 
     Args:
         f: A vector field, e.g. `[f_x, f_y, f_z]` where each `f_` component is a 1- to
-           3-D ndarray (`f_*` should all be the same size). Doesn't fail with `f=None`.
+            3-D ndarray (`f_*` should all be the same size). Doesn't fail with `f=None`.
 
     Returns:
         1D ndarray containing the linearized field (or `None`)
@@ -123,4 +123,3 @@ def unvec(
     if v is None:
         return None
     return v.reshape((nvdim, *shape), order='C')  # type: ignore
-

meanas/test/test_fdtd_pml.py

@@ -1,7 +1,10 @@
 import numpy
 import pytest
 
+from .. import fdtd
+from ..fdtd.base import maxwell_e, maxwell_h
 from ..fdtd.pml import cpml_params, updates_with_cpml
+from .utils import assert_close
 
 
 @pytest.mark.parametrize(
@@ -42,3 +45,202 @@ def test_updates_with_cpml_keeps_zero_fields_zero() -> None:
 
     assert not e.any()
     assert not h.any()
+
+
+def _unit_dxes(shape: tuple[int, int, int, int]) -> list[list[numpy.ndarray]]:
+    return [[numpy.ones(n, dtype=float) for n in shape[1:]] for _ in range(2)]
+
+
+def _real_field(shape: tuple[int, int, int, int], start: float) -> numpy.ndarray:
+    total = numpy.prod(shape, dtype=int)
+    return numpy.arange(start, start + total, dtype=float).reshape(shape) / total
+
+
+def _complex_field(shape: tuple[int, int, int, int], start: float) -> numpy.ndarray:
+    real = _real_field(shape, start)
+    imag = _real_field(shape, start + real.size)
+    return real + 1j * imag
+
+
+def test_updates_with_cpml_matches_base_updates_when_all_faces_disabled() -> None:
+    shape = (3, 4, 5, 6)
+    epsilon = _real_field(shape, 1.0) + 2.0
+    mu = _real_field(shape, 4.0) + 1.5
+    e = _real_field(shape, 10.0)
+    h = _real_field(shape, 100.0)
+    dxes = _unit_dxes(shape)
+    params = [[None, None] for _ in range(3)]
+
+    update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
+    update_e_base = maxwell_e(dt=0.1, dxes=dxes)
+    update_h_base = maxwell_h(dt=0.1, dxes=dxes)
+
+    e_cpml = e.copy()
+    h_cpml = h.copy()
+    e_base = e.copy()
+    h_base = h.copy()
+
+    update_e_cpml(e_cpml, h_cpml, epsilon)
+    update_e_base(e_base, h_base, epsilon)
+    update_h_cpml(e_cpml, h_cpml, mu)
+    update_h_base(e_base, h_base, mu)
+
+    assert_close(e_cpml, e_base)
+    assert_close(h_cpml, h_base)
+
+
+def test_updates_with_cpml_matches_base_updates_with_complex_dtype_when_all_faces_disabled() -> None:
+    shape = (3, 4, 5, 6)
+    epsilon = _real_field(shape, 1.0) + 2.0
+    mu = _real_field(shape, 4.0) + 1.5
+    e = _complex_field(shape, 10.0)
+    h = _complex_field(shape, 100.0)
+    dxes = _unit_dxes(shape)
+    params = [[None, None] for _ in range(3)]
+
+    update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon, dtype=complex)
+    update_e_base = maxwell_e(dt=0.1, dxes=dxes)
+    update_h_base = maxwell_h(dt=0.1, dxes=dxes)
+
+    e_cpml = e.copy()
+    h_cpml = h.copy()
+    e_base = e.copy()
+    h_base = h.copy()
+
+    update_e_cpml(e_cpml, h_cpml, epsilon)
+    update_e_base(e_base, h_base, epsilon)
+    update_h_cpml(e_cpml, h_cpml, mu)
+    update_h_base(e_base, h_base, mu)
+
+    assert_close(e_cpml, e_base)
+    assert_close(h_cpml, h_base)
+
+
+def test_updates_with_cpml_only_changes_the_configured_face_region() -> None:
+    shape = (3, 6, 6, 6)
+    epsilon = numpy.ones(shape, dtype=float)
+    mu = numpy.ones(shape, dtype=float)
+    e = _real_field(shape, 1.0)
+    h = _real_field(shape, 100.0)
+    dxes = _unit_dxes(shape)
+    thickness = 3
+
+    params = [[None, None] for _ in range(3)]
+    params[0][0] = cpml_params(axis=0, polarity=-1, dt=0.1, thickness=thickness)
+
+    update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
+    update_e_base = maxwell_e(dt=0.1, dxes=dxes)
+    update_h_base = maxwell_h(dt=0.1, dxes=dxes)
+
+    e_cpml = e.copy()
+    h_cpml = h.copy()
+    e_base = e.copy()
+    h_base = h.copy()
+
+    update_e_cpml(e_cpml, h_cpml, epsilon)
+    update_e_base(e_base, h_base, epsilon)
+    update_h_cpml(e_cpml, h_cpml, mu)
+    update_h_base(e_base, h_base, mu)
+
+    e_untouched = slice(thickness, None)
+    h_untouched = slice(thickness, -1)
+    assert_close(e_cpml[:, e_untouched, :, :], e_base[:, e_untouched, :, :])
+    assert_close(h_cpml[:, h_untouched, :, :], h_base[:, h_untouched, :, :])
+
+    changed_e = numpy.any(numpy.abs(e_cpml[:, :thickness, :, :] - e_base[:, :thickness, :, :]) > 1e-12)
+    changed_h = numpy.any(numpy.abs(h_cpml[:, :thickness, :, :] - h_base[:, :thickness, :, :]) > 1e-12)
+    assert changed_e
+    assert changed_h
+
+
+def test_cpml_plane_wave_phasor_decays_monotonically_through_outgoing_pml() -> None:
+    dt = 0.4
+    period_steps = 24
+    omega = 2 * numpy.pi / (period_steps * dt)
+    shape = (3, 80, 1, 1)
+    thickness = 8
+    source_x = 16
+    warmup_periods = 10
+    accumulation_periods = 6
+    total_steps = period_steps * (warmup_periods + accumulation_periods)
+
+    epsilon = numpy.ones(shape, dtype=float)
+    dxes = _unit_dxes(shape)
+    params = [[None, None] for _ in range(3)]
+    for polarity_index, polarity in enumerate((-1, 1)):
+        params[0][polarity_index] = cpml_params(axis=0, polarity=polarity, dt=dt, thickness=thickness)
+
+    update_e, update_h = updates_with_cpml(params, dt=dt, dxes=dxes, epsilon=epsilon)
+
+    e = numpy.zeros(shape, dtype=float)
+    h = numpy.zeros(shape, dtype=float)
+    e_accumulator = numpy.zeros((1, *shape), dtype=complex)
+
+    for step in range(total_steps):
+        update_e(e, h, epsilon)
+
+        source = numpy.cos(omega * (step + 0.5) * dt)
+        e[1, source_x, 0, 0] -= dt * source
+
+        if step >= period_steps * warmup_periods:
+            fdtd.accumulate_phasor_e(e_accumulator, omega, dt, e, step + 1)
+
+        update_h(e, h)
+
+    profile = numpy.abs(e_accumulator[0, 1, :, 0, 0])
+    right_pml = profile[-thickness:]
+    interior = profile[-thickness - 6:-thickness]
+    interior_level = interior.mean()
+
+    assert interior_level > 1.0
+    assert right_pml[-1] < interior_level / 100
+    assert profile[0] < interior_level / 100
+    assert numpy.all(numpy.diff(right_pml) <= interior_level * 1e-3)
+
+
+def test_cpml_point_source_total_energy_reaches_late_time_plateau() -> None:
+    dt = 0.3
+    period_steps = 24
+    omega = 2 * numpy.pi / (period_steps * dt)
+    cycles = 1000
+    sample_every_cycles = 50
+    sample_stride = period_steps * sample_every_cycles
+    shape = (3, 9, 9, 9)
+    thickness = 3
+    center = shape[1] // 2
+
+    epsilon = numpy.ones(shape, dtype=float)
+    dxes = _unit_dxes(shape)
+    params = [[None, None] for _ in range(3)]
+    for axis in range(3):
+        for polarity_index, polarity in enumerate((-1, 1)):
+            params[axis][polarity_index] = cpml_params(axis=axis, polarity=polarity, dt=dt, thickness=thickness)
+
+    update_e, update_h = updates_with_cpml(params, dt=dt, dxes=dxes, epsilon=epsilon)
+
+    e = numpy.zeros(shape, dtype=float)
+    h = numpy.zeros(shape, dtype=float)
+    sampled_energies: list[float] = []
+
+    for step in range(period_steps * cycles):
+        h_before = h.copy()
+        update_e(e, h, epsilon)
+
+        source = numpy.cos(omega * (step + 0.5) * dt)
+        e[1, center, center, center] -= dt * source
+
+        update_h(e, h)
+
+        if (step + 1) % sample_stride == 0:
+            total_energy = fdtd.energy_estep(h0=h_before, e1=e, h2=h, epsilon=epsilon, dxes=dxes).sum().real
+            sampled_energies.append(total_energy)
+
+    energies = numpy.asarray(sampled_energies)
+    late_window = energies[-5:]
+    previous_window = energies[-10:-5]
+    late_mean = late_window.mean()
+
+    assert energies.size == cycles // sample_every_cycles
+    assert late_mean > 0.1
+    assert (late_window.max() - late_window.min()) / late_mean < 1e-4
+    assert abs(late_mean - previous_window.mean()) / late_mean < 1e-4

mkdocs.yml

@@ -10,6 +10,19 @@ strict: false
 theme:
   name: material
   font: false
+  palette:
+    - scheme: slate
+      primary: blue grey
+      accent: cyan
+      toggle:
+        icon: material/weather-sunny
+        name: Switch to light mode
+    - scheme: default
+      primary: teal
+      accent: indigo
+      toggle:
+        icon: material/weather-night
+        name: Switch to dark mode
   features:
     - navigation.indexes
     - navigation.sections

pyproject.toml

@@ -63,9 +63,11 @@ docs = [
     "mkdocs-print-site-plugin>=2.3",
     "pymdown-extensions>=10.7",
     "htmlark>=1.0",
+    "ruff>=0.6",
     ]
 examples = [
     "matplotlib>=3.10.8",
+    "scikit-rf>=1.0",
 ]
 test = ["pytest", "coverage"]