diff --git a/README.md b/README.md index b400796..ce9bcc1 100644 --- a/README.md +++ b/README.md @@ -163,11 +163,6 @@ 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. diff --git a/docs/index.md b/docs/index.md index 5475af8..bc4cdeb 100644 --- a/docs/index.md +++ b/docs/index.md @@ -24,10 +24,6 @@ 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 diff --git a/docs/stylesheets/extra.css b/docs/stylesheets/extra.css index 090b5e7..d70a077 100644 --- a/docs/stylesheets/extra.css +++ b/docs/stylesheets/extra.css @@ -11,42 +11,3 @@ .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; -} diff --git a/examples/eme.py b/examples/eme.py deleted file mode 100644 index 6215cbc..0000000 --- a/examples/eme.py +++ /dev/null @@ -1,217 +0,0 @@ -""" -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() diff --git a/examples/eme_bend.py b/examples/eme_bend.py deleted file mode 100644 index caff4df..0000000 --- a/examples/eme_bend.py +++ /dev/null @@ -1,310 +0,0 @@ -""" -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() diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py index 6425a29..18aaade 100644 --- a/meanas/fdfd/operators.py +++ b/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). - The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + as containing a perfect electrical conductor (PEC). + 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). - The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + as containing a perfect magnetic conductor (PMC). + 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). - The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + as containing a perfect electrical conductor (PEC). + 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). - The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + as containing a perfect magnetic conductor (PMC). + 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). - The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + as containing a perfect electrical conductor (PEC). + 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). - The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + as containing a perfect magnetic conductor (PMC). + 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). - The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + as containing a perfect magnetic conductor (PMC). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) Returns: Sparse matrix for converting E to H. @@ -483,3 +483,4 @@ def e_boundary_source( # (numpy.roll(mask, +1, axis=2) != mask)) return sparse.diags_array(jmask.astype(int)) @ full + diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py index e7dfd22..66cc7cc 100644 --- a/meanas/fdfd/waveguide_3d.py +++ b/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: diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index 201f709..caedfaf 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -216,13 +216,13 @@ def solve_modes( of the bent waveguide with the specified mode number. Args: - mode_numbers: Mode numbers to solve, 0-indexed. + mode_number: Number of the mode, 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. diff --git a/meanas/fdmath/operators.py b/meanas/fdmath/operators.py index 0c64ae7..19ccb80 100644 --- a/meanas/fdmath/operators.py +++ b/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. diff --git a/meanas/fdmath/vectorization.py b/meanas/fdmath/vectorization.py index 77b722f..44d8b74 100644 --- a/meanas/fdmath/vectorization.py +++ b/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,3 +123,4 @@ def unvec( if v is None: return None return v.reshape((nvdim, *shape), order='C') # type: ignore + diff --git a/meanas/test/test_fdtd_pml.py b/meanas/test/test_fdtd_pml.py index 9d8aef8..6d118c6 100644 --- a/meanas/test/test_fdtd_pml.py +++ b/meanas/test/test_fdtd_pml.py @@ -1,10 +1,7 @@ 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( @@ -45,202 +42,3 @@ 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 diff --git a/mkdocs.yml b/mkdocs.yml index 379af9b..a0dcd2c 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -10,19 +10,6 @@ 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 diff --git a/pyproject.toml b/pyproject.toml index 013631a..4b81fa6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -63,11 +63,9 @@ 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"]