From 7e407dd1ac6ecb5b8b4484c8d6a8e26a7c5b5f77 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Sun, 18 Oct 2026 23:07:53 -0700 Subject: [PATCH] [docs] document scattering conventions and frozen-basis workflows Update examples for sampled-interface tapers, explicit physical CPML spacing, copied SCPML widths, and the consolidated complex mode solver. Document breaking API changes, overlap measures, and PML convergence. --- README.md | 15 +++--- docs/api/fdfd.md | 14 +++++ docs/api/fdtd.md | 5 ++ docs/api/waveguides.md | 86 +++++++++++++++++++++++++++++- docs/index.md | 2 +- examples/eme_taper_cmt.py | 11 ++-- examples/eme_taper_pml_modes.py | 38 +++++++++---- examples/fdtd.py | 4 +- examples/waveguide.py | 4 +- meanas/test/test_examples_smoke.py | 2 +- 10 files changed, 152 insertions(+), 29 deletions(-) diff --git a/README.md b/README.md index 0061688..d1b7f95 100644 --- a/README.md +++ b/README.md @@ -49,7 +49,7 @@ linear systems, ideally with double precision. **Requirements:** -* python >=3.11 +* python >=3.12 * numpy * scipy @@ -118,11 +118,12 @@ python3 -m pytest -q | tee test_results.txt - `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.fdfd.waveguide_2d_complex`: targeted signed-`beta` waveguide modes, - transpose-left / Lorentz-dual projection fields, PML-mode diagnostics, and - overlap-based tracking for open cross-sections. + solvers, including targeted complex bases from the default half-sized + beta-squared or explicit signed-beta formulations, transpose-left and + Lorentz-dual projection fields on a shared field-array axis, separately + returned PML-mode diagnostics, overlap-based tracking, + 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, @@ -175,7 +176,7 @@ The tracked examples under `examples/` are the intended entry points for users: - `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/eme_taper_cmt.py`: sampled cross-section local-mode CMT for a +- `examples/eme_taper_cmt.py`: sampled-interface EME with stable S cascading for a continuously varying rib-waveguide taper. - `examples/eme_taper_pml_modes.py`: targeted complex/PML mode solving, tracking, and construction of EME taper sections. diff --git a/docs/api/fdfd.md b/docs/api/fdfd.md index e47c278..83328f0 100644 --- a/docs/api/fdfd.md +++ b/docs/api/fdfd.md @@ -13,3 +13,17 @@ ::: meanas.fdfd.scpml ::: meanas.fdfd.farfield + +## SCPML conventions + +`stretch_with_scpml` now uses `-1` for the low-index face and `+1` for the +high-index face, matching CPML and propagation conventions. It returns copied +widths by default; pass `copy=False` to request mutation explicitly. + +Both SCPML builders accept `kappa_function(depth)` for real stretching and +`max_imag_ratio` to limit the final dimensionless `abs(Im(s))/Re(s)`. Start a +mode-solver convergence study with `max_imag_ratio=1.0`: overly strong imaginary +stretching can produce above-light-line PML artifacts. The cap is not a general +spectral guarantee, and weaker absorption can increase boundary reflections. +Check domain size, PML thickness, strength, mode localization, and residuals. +The general scattering profile remains uncapped by default. diff --git a/docs/api/fdtd.md b/docs/api/fdtd.md index 3e1823c..f3e7ac2 100644 --- a/docs/api/fdtd.md +++ b/docs/api/fdtd.md @@ -13,3 +13,8 @@ ::: meanas.fdtd.energy ::: meanas.fdtd.phasor + +CPML derivatives use electric widths for forward differences and magnetic +widths for backward differences. `cpml_params(..., dx=1.0)` specifies the uniform +physical cell width normal to a face; arbitrary nonuniform PML grading is not +implemented by this constructor. Scale the timestep consistently with the grid. diff --git a/docs/api/waveguides.md b/docs/api/waveguides.md index a97e618..464baba 100644 --- a/docs/api/waveguides.md +++ b/docs/api/waveguides.md @@ -2,8 +2,90 @@ ::: meanas.fdfd.waveguide_2d -::: meanas.fdfd.waveguide_2d_complex - ::: meanas.fdfd.waveguide_3d ::: meanas.fdfd.waveguide_cyl + +## Overlap measures and normalization + +`waveguide_2d.inner_product(..., measure="physical", conj_h=True)` uses real +Yee widths for physical Poynting power. `measure="stretched"` keeps complex +widths for coordinate reaction. `lorentz_metric` and `lorentz_overlap` default +to the stretched measure. Use the same grid and measure when constructing and +using dual fields. A real overlap grid still defines a consistent constructed +dual basis, even when the eigenproblem uses a separate PML grid. + +`normalized_fields_e/h(..., normalization="lorentz")` normalizes unconjugated +self-reaction to one. This fixes phase up to a sign and is distinct from unit +forward power, which remains the default. Zero-power or self-orthogonal modes +raise `ValueError`; open self-orthogonal modes require paired primal/dual +normalization through `ModeBasis`. + +## Targeted and reproducible mode solving + +For modes near an effective index guess, use `solve_modes_near` with +`target_wavenumbers=[omega * n_guess]`. It targets the electric beta-squared +operator at `(omega * n_guess)**2`; `solve_modes_near_beta` instead targets the +first-order signed-beta operator. Both use the LU shift-invert backend and +support `eigensolver_opts={"v0": starting_vector}`. Factorization cost and +memory depend on the grid, so this is not universally faster than indexed +solving. `solve_modes` keeps its global mode-number semantics and now accepts +`rng` or `v0` for reproducible initialization. + +`waveguide_3d.mode_from_2d` accepts a supplied transverse electric eigenvector +and reuses field reconstruction, normalization, numerical-dispersion +correction, Yee staggering, rotation, and full-grid embedding without solving +an eigenproblem. + +## Interface and sampled-taper scattering + +All scattering matrices use outgoing rows and incoming columns, ordered +`[left, right]`: + +```text +[b_left ] [R_left T_right_to_left] [a_left ] +[b_right] = [T_left_to_right R_right] [a_right] +``` + +`get_s` and `get_abcd` require explicit `dxes`. Mixed Petrov-Galerkin testing is +the default: left electric duals test H continuity, and right magnetic duals +test E continuity. Both incidence directions are solved together with full +projection matrices. Mixed scattering supports unequal mode counts. +`formulation="galerkin"` selects equal-count left testing of both equations. +Neither formulation silently clips transmission. The optional `force_nogain` +and `force_reciprocal` flags are numerical projections, not general physical +tests for arbitrary open or evanescent modal coordinates. + +`get_taper_s` replaces finite-difference CMT with half propagation in the left +basis, a midpoint interface, and half propagation in the right basis for each +sampled interval. Redheffer composition avoids growing passive backward +transfer branches. Refine section sampling and modal/PML bases to establish +convergence. The ineffective `rtol`, `atol`, and `max_step` arguments have been +removed. The existing `eme_taper_cmt.py` example filename is retained. + +```python +prepared = eme.prepare_taper(sections, dxes) +s_original = eme.evaluate_taper_s(prepared) +s_longer = eme.evaluate_taper_s(prepared, z_coords=prepared.z_coords * 1.2) +s_changed_beta = eme.evaluate_taper_s( + prepared, wavenumber_matrices=beta_matrices, +) +``` + +`PreparedTaper` owns read-only copies of phase-aligned interface data. Beta +matrix overrides have shape `(section, mode, mode)` in the original supplied +section bases; recorded basis transforms are applied automatically. Explicit +backward branches require backward beta overrides when forward beta is changed. +Evaluation performs no field overlaps or eigensolves. Length sweeps reuse the +sampled model; wavelength sweeps with reused fields are explicitly frozen-basis +approximations. Reprepare when fields, duals, grids, or formulations change. + +`cascade_s(left, right, internal_modes=...)` connects spatially ordered S +matrices, including unequal external ports. Public `abcd_to_s` and `s_to_abcd` +use `[f_right, b_right] = ABCD @ [f_left, b_left]` and require equal port counts. +`get_taper_abcd` converts the stable S result; zero transmission or an +unrepresentable transfer matrix raises `LinAlgError`. Prefer S for long lossy +devices. Singular interface and cascade systems also raise contextual errors +rather than silently substituting pseudoinverses. + +::: meanas.fdfd.eme diff --git a/docs/index.md b/docs/index.md index 1b5ca6b..22c9379 100644 --- a/docs/index.md +++ b/docs/index.md @@ -28,7 +28,7 @@ Relevant starting examples: 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/eme_taper_cmt.py` for local-mode CMT through sampled continuously +- `examples/eme_taper_cmt.py` for sampled-interface EME through continuously varying taper cross-sections - `examples/fdfd.py` for direct frequency-domain waveguide excitation diff --git a/examples/eme_taper_cmt.py b/examples/eme_taper_cmt.py index c5a8f1c..c28e68b 100644 --- a/examples/eme_taper_cmt.py +++ b/examples/eme_taper_cmt.py @@ -1,5 +1,5 @@ """ -Local-mode CMT example for a continuously varying rib-waveguide taper. +Sampled-interface EME example for a continuously varying rib-waveguide taper. This example keeps geometry construction outside `meanas.fdfd.eme`: it samples a width taper at several cross-sections, solves and normalizes each local mode with @@ -112,14 +112,15 @@ def print_summary( print('sampled taper effective indices:', ', '.join(f'{value:.5f}' for value in neffs)) print(f'abrupt endpoint reflection |S_00|^2 = {abs(abrupt_s[0, 0]) ** 2:.6f}') print(f'abrupt endpoint transmission |S_{n_modes},0|^2 = {abs(abrupt_s[n_modes, 0]) ** 2:.6f}') - print(f'taper CMT reflection |S_00|^2 = {abs(taper_s[0, 0]) ** 2:.6f}') - print(f'taper CMT transmission |S_{n_modes},0|^2 = {abs(taper_s[n_modes, 0]) ** 2:.6f}') - print(f'taper CMT total output power = {numpy.sum(abs(taper_s[:, 0]) ** 2):.6f}') + print(f'taper reflection |S_00|^2 = {abs(taper_s[0, 0]) ** 2:.6f}') + print(f'taper transmission |S_{n_modes},0|^2 = {abs(taper_s[n_modes, 0]) ** 2:.6f}') + print(f'taper total output power = {numpy.sum(abs(taper_s[:, 0]) ** 2):.6f}') def main() -> None: sections, neffs, _omega, dxes_2d = solve_taper_sections() - taper_s = eme.get_taper_s(sections, dxes=dxes_2d) + prepared = eme.prepare_taper(sections, dxes=dxes_2d) + taper_s = eme.evaluate_taper_s(prepared) abrupt_s = eme.get_s( sections[0].modes, sections[0].wavenumbers, diff --git a/examples/eme_taper_pml_modes.py b/examples/eme_taper_pml_modes.py index cb71d30..9e37483 100644 --- a/examples/eme_taper_pml_modes.py +++ b/examples/eme_taper_pml_modes.py @@ -5,7 +5,7 @@ from __future__ import annotations import numpy from meanas.fdmath import vec -from meanas.fdfd import eme, scpml, waveguide_2d_complex +from meanas.fdfd import eme, scpml, waveguide_2d WAVELENGTH = 1.55 @@ -34,6 +34,7 @@ def build_grids() -> tuple[list[list[numpy.ndarray]], list[list[numpy.ndarray]]] omega=OMEGA, epsilon_effective=EPS_CLADDING, thickness=PML_CELLS, + max_imag_ratio=1.0, ) return physical, stretched @@ -61,29 +62,48 @@ def main() -> None: """Solve, track, and summarize the sampled complex modal bases.""" physical_dxes, solve_dxes = build_grids() z_coords = numpy.array([0.0, 5.0, 10.0]) - bases = [ - waveguide_2d_complex.solve_modes_near( + epsilon_samples = [build_epsilon(edge_fill) for edge_fill in (0.30, 0.35, 0.40)] + solved = [ + waveguide_2d.solve_modes_near( OMEGA, solve_dxes, - build_epsilon(edge_fill), + epsilon, target_wavenumbers=[1.45 * OMEGA], modes_per_target=3, candidate_margin=3, overlap_dxes=physical_dxes, ) - for edge_fill in (0.30, 0.35, 0.40) + for epsilon in epsilon_samples + ] + bases = [basis for basis, _diagnostics in solved] + tracked = waveguide_2d.track_mode_path(bases, overlap_dxes=physical_dxes) + tracked_diagnostics = [ + waveguide_2d.compute_mode_diagnostics( + basis, + OMEGA, + solve_dxes, + epsilon, + target_wavenumbers=[1.45 * OMEGA], + overlap_dxes=physical_dxes, + ) + for basis, epsilon in zip(tracked, epsilon_samples, strict=True) ] - tracked = waveguide_2d_complex.track_mode_path(bases) sections = eme.mode_sections_from_bases(z_coords, tracked) + prepared = eme.prepare_taper(sections, physical_dxes) + # Fields and interface overlaps are reused: these are frozen-basis length + # sweeps. Wavelength-dependent beta overrides would also freeze the fields. + length_sweep = [eme.evaluate_taper_s(prepared, z_coords=z_coords * scale) + for scale in (0.8, 1.0, 1.2)] - for z_coord, basis in zip(z_coords, tracked, strict=True): + for z_coord, basis, diagnostics in zip(z_coords, tracked, tracked_diagnostics, strict=True): print( f'z={z_coord:5.1f}: ' f'beta={basis.wavenumbers[0]:.6g}, ' - f'PML participation={basis.diagnostics.pml_participation[0]:.3f}, ' - f'right residual={basis.diagnostics.right_residuals[0]:.2e}', + f'PML participation={diagnostics.pml_participation[0]:.3f}, ' + f'right residual={diagnostics.right_residuals[0]:.2e}', ) print(f'constructed {len(sections)} tracked EME sections with {len(sections[0].modes)} modes each') + print('frozen-basis length sweep |T_00|:', [float(abs(ss[len(tracked[0]), 0])) for ss in length_sweep]) if __name__ == '__main__': diff --git a/examples/fdtd.py b/examples/fdtd.py index fd6026d..de7b33d 100644 --- a/examples/fdtd.py +++ b/examples/fdtd.py @@ -145,7 +145,7 @@ def main() -> None: # PMLs in every direction pml_params = [ - [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_air ** 2) + [cpml_params(axis=dd, polarity=pp, dt=dt, dx=dx, thickness=pml_thickness, epsilon_eff=n_air ** 2) for pp in (-1, +1)] for dd in range(3)] update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, dxes=dxes, epsilon=epsilon, dtype=complex) @@ -203,7 +203,7 @@ def main() -> None: dxes_fdfd = copy.deepcopy(dxes) 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) + dxes_fdfd = 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) diff --git a/examples/waveguide.py b/examples/waveguide.py index f243924..7301c6c 100644 --- a/examples/waveguide.py +++ b/examples/waveguide.py @@ -273,7 +273,7 @@ def main2(): # PMLs in every direction pml_params = [ - [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_cladding ** 2) + [cpml_params(axis=dd, polarity=pp, dt=dt, dx=dx, thickness=pml_thickness, epsilon_eff=n_cladding ** 2) for pp in (-1, +1)] for dd in range(3)] update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, dxes=dxes, epsilon=epsilon, dtype=complex) @@ -330,7 +330,7 @@ def main2(): dxes_fdfd = copy.deepcopy(dxes) 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) + dxes_fdfd = 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) diff --git a/meanas/test/test_examples_smoke.py b/meanas/test/test_examples_smoke.py index feefcb3..b418390 100644 --- a/meanas/test/test_examples_smoke.py +++ b/meanas/test/test_examples_smoke.py @@ -57,4 +57,4 @@ def test_eme_taper_cmt_example_smoke_runs(tmp_path: Path) -> None: assert result.returncode == 0, result.stdout + result.stderr assert 'sampled taper effective indices:' in result.stdout - assert 'taper CMT transmission' in result.stdout + assert 'taper transmission' in result.stdout