From 0e89fe0460480770c0cb0197490a83b10c54f603 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Tue, 14 Jul 2026 21:15:45 -0700 Subject: [PATCH] [fdfd.waveguide_2d_complex] add 2d complex modesolver --- README.md | 5 + docs/api/waveguides.md | 2 + examples/eme_taper_pml_modes.py | 90 +++ meanas/eigensolvers.py | 135 +++- meanas/fdfd/__init__.py | 1 + meanas/fdfd/eme.py | 200 ++++- meanas/fdfd/waveguide_2d_complex.py | 881 ++++++++++++++++++++++ meanas/test/test_eigensolvers_numerics.py | 27 +- meanas/test/test_waveguide_2d_complex.py | 694 +++++++++++++++++ 9 files changed, 2021 insertions(+), 14 deletions(-) create mode 100644 examples/eme_taper_pml_modes.py create mode 100644 meanas/fdfd/waveguide_2d_complex.py create mode 100644 meanas/test/test_waveguide_2d_complex.py diff --git a/README.md b/README.md index d179c9d..0061688 100644 --- a/README.md +++ b/README.md @@ -120,6 +120,9 @@ python3 -m pytest -q | tee test_results.txt - `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. - `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, @@ -174,6 +177,8 @@ The tracked examples under `examples/` are the intended entry points for users: example built with `scikit-rf`. - `examples/eme_taper_cmt.py`: sampled cross-section local-mode CMT 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. - `examples/fdfd.py`: direct frequency-domain waveguide excitation and overlap / Poynting analysis without a time-domain run. diff --git a/docs/api/waveguides.md b/docs/api/waveguides.md index 834614c..a97e618 100644 --- a/docs/api/waveguides.md +++ b/docs/api/waveguides.md @@ -2,6 +2,8 @@ ::: meanas.fdfd.waveguide_2d +::: meanas.fdfd.waveguide_2d_complex + ::: meanas.fdfd.waveguide_3d ::: meanas.fdfd.waveguide_cyl diff --git a/examples/eme_taper_pml_modes.py b/examples/eme_taper_pml_modes.py new file mode 100644 index 0000000..cb71d30 --- /dev/null +++ b/examples/eme_taper_pml_modes.py @@ -0,0 +1,90 @@ +"""Build tracked complex/PML modal bases for a slowly varying taper.""" + +from __future__ import annotations + +import numpy + +from meanas.fdmath import vec +from meanas.fdfd import eme, scpml, waveguide_2d_complex + + +WAVELENGTH = 1.55 +OMEGA = 2 * numpy.pi / WAVELENGTH +SHAPE = (17, 17) +SPACING = 0.25 +PML_CELLS = 3 +EPS_CLADDING = 1.45 ** 2 +EPS_CORE = 2.0 ** 2 + + +def build_grids() -> tuple[list[list[numpy.ndarray]], list[list[numpy.ndarray]]]: + """Build matching physical and complex-stretched transverse Yee grids. + + Returns: + The real overlap grid followed by the PML-stretched solve grid. + """ + physical = [[numpy.full(size, SPACING) for size in SHAPE] for _ in range(2)] + stretched = [[values.copy() for values in family] for family in physical] + for axis in (0, 1): + for polarity in (-1, +1): + stretched = scpml.stretch_with_scpml( + stretched, + axis=axis, + polarity=polarity, + omega=OMEGA, + epsilon_effective=EPS_CLADDING, + thickness=PML_CELLS, + ) + return physical, stretched + + +def build_epsilon(edge_fill: float) -> numpy.ndarray: + """Construct and vectorize one taper cross-section permittivity profile. + + Args: + edge_fill: Fractional core fill used for the two cells beside the + central three-cell-wide core. + + Returns: + Vectorized three-component relative permittivity field. + """ + epsilon = numpy.full((3, *SHAPE), EPS_CLADDING) + center_x, center_y = (size // 2 for size in SHAPE) + epsilon[:, center_x - 1:center_x + 2, center_y] = EPS_CORE + mixed_epsilon = EPS_CLADDING + edge_fill * (EPS_CORE - EPS_CLADDING) + epsilon[:, center_x - 1:center_x + 2, center_y - 1] = mixed_epsilon + epsilon[:, center_x - 1:center_x + 2, center_y + 1] = mixed_epsilon + return vec(epsilon) + + +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( + OMEGA, + solve_dxes, + build_epsilon(edge_fill), + 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) + ] + tracked = waveguide_2d_complex.track_mode_path(bases) + sections = eme.mode_sections_from_bases(z_coords, tracked) + + for z_coord, basis in zip(z_coords, tracked, 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}', + ) + print(f'constructed {len(sections)} tracked EME sections with {len(sections[0].modes)} modes each') + + +if __name__ == '__main__': + main() diff --git a/meanas/eigensolvers.py b/meanas/eigensolvers.py index 98e7a15..b00894e 100644 --- a/meanas/eigensolvers.py +++ b/meanas/eigensolvers.py @@ -2,11 +2,142 @@ Solvers for eigenvalue / eigenvector problems """ from collections.abc import Callable +import dataclasses import numpy from numpy.typing import NDArray, ArrayLike from numpy.linalg import norm from scipy import sparse import scipy.sparse.linalg as spalg +from scipy.optimize import linear_sum_assignment + + +@dataclasses.dataclass(frozen=True) +class TargetedEigensolution: + """Paired right and transpose-left eigenvectors near a spectral target. + + Attributes: + eigenvalues: Eigenvalues associated with the columns of both vector + matrices. + right_eigenvectors: Right eigenvectors stored column-wise. + left_eigenvectors: Bilinear-left eigenvectors stored column-wise. These + satisfy a transpose, rather than Hermitian-transpose, eigenproblem. + """ + + eigenvalues: NDArray[numpy.complex128] + right_eigenvectors: NDArray[numpy.complex128] + left_eigenvectors: NDArray[numpy.complex128] + + +def scipy_eigs_near( + operator: sparse.sparray, + target: complex, + how_many: int, + *, + tolerance: float = 1e-10, + max_iterations: int | None = None, + ncv: int | None = None, + pairing_rtol: float = 1e-7, + v0: ArrayLike | None = None, + ) -> TargetedEigensolution: + """ + Solve for paired right and bilinear-left eigenvectors near ``target``. + + One sparse LU factorization of ``operator - target * I`` is reused for the + right solve and for the transpose-left solve. The left vectors therefore + satisfy ``left.T @ operator = eigenvalue * left.T``; this deliberately uses + a transpose rather than a Hermitian transpose so the result is suitable for + unconjugated Lorentz-reciprocity projections. + + Args: + operator: Square sparse operator to diagonalize. + target: Shift-invert spectral target. + how_many: Number of eigenpairs to return. + tolerance: ARPACK convergence tolerance for both solves. + max_iterations: Optional ARPACK iteration limit. + ncv: Optional number of Arnoldi vectors. A size-aware default is used + when omitted. + pairing_rtol: Relative tolerance for matching right and transpose-left + eigenvalues. + v0: Optional finite starting vector shared by both solves. + + Returns: + Paired eigenvalues and column-wise right and transpose-left vectors. + + Raises: + ValueError: If the dimensions, counts, tolerances, or starting vector + are invalid. + RuntimeError: If the independently solved spectra cannot be paired. + """ + if how_many < 1: + raise ValueError('how_many must be positive') + if not numpy.isfinite(tolerance) or tolerance <= 0: + raise ValueError('tolerance must be finite and positive') + if not numpy.isfinite(pairing_rtol) or pairing_rtol <= 0: + raise ValueError('pairing_rtol must be finite and positive') + + operator = sparse.csc_array(operator, dtype=complex) + size = operator.shape[0] + if operator.shape != (size, size): + raise ValueError('operator must be square') + if how_many >= size - 1: + raise ValueError('how_many must be smaller than operator dimension minus one') + + shifted = operator - target * sparse.eye_array(size, format='csc', dtype=complex) + factor = spalg.splu(shifted) + dtype = numpy.dtype(numpy.complex128) + inverse = spalg.LinearOperator( + operator.shape, + matvec=lambda vector: factor.solve(vector, trans='N'), + dtype=dtype, + ) + inverse_transpose = spalg.LinearOperator( + operator.shape, + matvec=lambda vector: factor.solve(vector, trans='T'), + dtype=dtype, + ) + if ncv is None: + ncv = min(size, max(2 * how_many + 1, 4 * how_many)) + start = None if v0 is None else numpy.asarray(v0, dtype=complex) + if start is not None and (start.shape != (size,) or not numpy.isfinite(start).all()): + raise ValueError('v0 must be a finite vector matching the operator dimension') + + eigenvalues, right = spalg.eigs( + operator, + k=how_many, + sigma=target, + which='LM', + OPinv=inverse, + tol=tolerance, + maxiter=max_iterations, + ncv=ncv, + v0=start, + ) + left_values, left = spalg.eigs( + operator.T, + k=how_many, + sigma=target, + which='LM', + OPinv=inverse_transpose, + tol=tolerance, + maxiter=max_iterations, + ncv=ncv, + v0=start, + ) + + rows, columns = linear_sum_assignment(numpy.abs(eigenvalues[:, None] - left_values[None, :])) + paired_columns = columns[numpy.argsort(rows)] + left_values = left_values[paired_columns] + left = left[:, paired_columns] + pairing_error = numpy.abs(eigenvalues - left_values) + pairing_scale = numpy.maximum(1.0, numpy.maximum(numpy.abs(eigenvalues), numpy.abs(left_values))) + if numpy.any(pairing_error > pairing_rtol * pairing_scale): + raise RuntimeError('could not pair right and transpose-left eigenvalues') + + return TargetedEigensolution( + eigenvalues=numpy.asarray(eigenvalues, dtype=complex), + right_eigenvectors=numpy.asarray(right, dtype=complex), + left_eigenvectors=numpy.asarray(left, dtype=complex), + ) def power_iteration( @@ -67,12 +198,14 @@ def rayleigh_quotient_iteration( (operator - sparse.eye_array(operator.shape[0])) def shift(eigval: float) -> sparse.sparray: + """Construct a sparse scalar identity shift.""" return eigval * sparse.eye_array(operator.shape[0]) if solver is None: solver = spalg.spsolve except TypeError: def shift(eigval: float) -> spalg.LinearOperator: + """Construct a matrix-free scalar identity shift.""" return spalg.LinearOperator( shape=operator.shape, dtype=operator.dtype, @@ -80,6 +213,7 @@ def rayleigh_quotient_iteration( ) if solver is None: def solver(A: spalg.LinearOperator, b: ArrayLike) -> NDArray[numpy.complex128]: + """Solve a matrix-free shifted system with BiCGSTAB.""" return spalg.bicgstab(A, b)[0] assert solver is not None @@ -140,4 +274,3 @@ def signed_eigensolve( k = eigenvalues.argsort() return eigenvalues[k], eigenvectors[:, k] - diff --git a/meanas/fdfd/__init__.py b/meanas/fdfd/__init__.py index e94adc3..23139d2 100644 --- a/meanas/fdfd/__init__.py +++ b/meanas/fdfd/__init__.py @@ -96,6 +96,7 @@ from . import ( functional as functional, scpml as scpml, waveguide_2d as waveguide_2d, + waveguide_2d_complex as waveguide_2d_complex, waveguide_3d as waveguide_3d, ) # from . import farfield, bloch TODO diff --git a/meanas/fdfd/eme.py b/meanas/fdfd/eme.py index a9de4dd..016b439 100644 --- a/meanas/fdfd/eme.py +++ b/meanas/fdfd/eme.py @@ -24,6 +24,7 @@ compose into larger workflows. from collections.abc import Sequence import dataclasses +import operator import numpy from numpy.typing import NDArray from scipy import linalg @@ -31,6 +32,7 @@ from scipy import sparse from ..fdmath import dx_lists2_t, vcfdfield2 from .waveguide_2d import inner_product +from .waveguide_2d_complex import ComplexModeBasis type wavenumber_seq = Sequence[complex] | NDArray[numpy.complexfloating] | NDArray[numpy.floating] type mode_seq = Sequence[Sequence[vcfdfield2]] @@ -54,6 +56,10 @@ class ModeSection: dual_backward_modes: Optional backward dual / adjoint projection modes. If omitted, they are synthesized from `dual_modes` when available, otherwise from `backward_modes`. + wavenumber_matrix: Optional dense beta operator in the supplied forward + basis. If omitted, ``diag(wavenumbers)`` is used. + backward_wavenumber_matrix: Optional dense beta operator for explicit + backward modes. Synthesized reciprocal branches use ``-wavenumber_matrix``. """ z: float @@ -63,6 +69,8 @@ class ModeSection: backward_wavenumbers: wavenumber_seq | None = None dual_modes: mode_seq | None = None dual_backward_modes: mode_seq | None = None + wavenumber_matrix: NDArray[numpy.complex128] | None = None + backward_wavenumber_matrix: NDArray[numpy.complex128] | None = None def _validate_port_modes( @@ -70,6 +78,7 @@ def _validate_port_modes( ehs: mode_seq, wavenumbers: wavenumber_seq, ) -> tuple[tuple[int, ...], tuple[int, ...]]: + """Validate a nonempty modal list and return its common E/H shapes.""" if len(ehs) != len(wavenumbers): raise ValueError(f'{name} mode list and wavenumber list must have the same length') if not ehs: @@ -102,6 +111,7 @@ def _validate_dual_modes( reference_shape: tuple[int, ...], wavenumbers: wavenumber_seq, ) -> mode_seq | None: + """Validate optional dual modes against a primal modal field shape.""" if dual_ehs is None: return None @@ -115,6 +125,7 @@ def _as_wavenumber_array( name: str, wavenumbers: wavenumber_seq, ) -> NDArray[numpy.complex128]: + """Convert finite propagation constants to a one-dimensional array.""" array = numpy.asarray(wavenumbers, dtype=complex) if array.ndim != 1: raise ValueError(f'{name} must be a one-dimensional sequence') @@ -123,9 +134,26 @@ def _as_wavenumber_array( return array +def _as_wavenumber_matrix( + name: str, + matrix: NDArray[numpy.complex128] | None, + wavenumbers: NDArray[numpy.complex128], + ) -> NDArray[numpy.complex128]: + """Return a validated propagation matrix or the diagonal spectrum.""" + if matrix is None: + return numpy.diag(wavenumbers) + array = numpy.asarray(matrix, dtype=complex) + if array.shape != (wavenumbers.size, wavenumbers.size): + raise ValueError(f'{name} must have one row and column per mode') + if not numpy.isfinite(array).all(): + raise ValueError(f'{name} must contain only finite values') + return array + + def _as_mode_arrays( ehs: mode_seq, ) -> list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]: + """Copy a modal sequence into complex NumPy array views.""" return [ (numpy.asarray(e_field, dtype=complex), numpy.asarray(h_field, dtype=complex)) for e_field, h_field in ehs @@ -137,6 +165,7 @@ def _lorentz_overlap( mode_b: tuple[vcfdfield2, vcfdfield2], dxes: dx_lists2_t, ) -> complex: + """Evaluate the symmetric unconjugated Lorentz overlap of two modes.""" e_a, h_a = mode_a e_b, h_b = mode_b return 0.5 * ( @@ -150,6 +179,7 @@ def _lorentz_derivative_overlap( derivative_b: tuple[vcfdfield2, vcfdfield2], dxes: dx_lists2_t, ) -> complex: + """Project a modal derivative through the symmetric Lorentz form.""" e_a, h_a = mode_a de_b, dh_b = derivative_b return 0.5 * ( @@ -163,8 +193,10 @@ def _phase_align_modes( current: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], dxes: dx_lists2_t, previous_dual: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]] | None = None, - ) -> list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]: + ) -> tuple[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], NDArray[numpy.complex128]]: + """Align current modal phases and return their diagonal basis transform.""" aligned = [] + scales = [] test_modes = previous if previous_dual is None else previous_dual for index, (previous_mode, current_mode, test_mode) in enumerate(zip(previous, current, test_modes, strict=True)): overlap = _lorentz_overlap(test_mode, current_mode, dxes) @@ -175,8 +207,10 @@ def _phase_align_modes( raise ValueError(f'cannot phase-track mode {index}: mode dual-overlap is zero') phase = (overlap / abs(overlap)) / (self_overlap / abs(self_overlap)) e_field, h_field = current_mode - aligned.append((e_field / phase, h_field / phase)) - return aligned + scale = 1 / phase + scales.append(scale) + aligned.append((e_field * scale, h_field * scale)) + return aligned, numpy.diag(scales) def _section_branches( @@ -192,14 +226,18 @@ def _section_branches( NDArray[numpy.complex128], list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], + NDArray[numpy.complex128], + NDArray[numpy.complex128], tuple[int, ...], ]: + """Validate one section and construct its forward and backward branches.""" z_coord = float(section.z) if not numpy.isfinite(z_coord): raise ValueError(f'sections[{index}].z must be finite') shape, _h_shape = _validate_port_modes(f'sections[{index}].modes', section.modes, section.wavenumbers) wavenumbers = _as_wavenumber_array(f'sections[{index}].wavenumbers', section.wavenumbers) + wavenumber_matrix = _as_wavenumber_matrix(f'sections[{index}].wavenumber_matrix', section.wavenumber_matrix, wavenumbers) if expected_count is not None and len(wavenumbers) != expected_count: raise ValueError('all taper sections must contain the same number of modes') if expected_shape is not None and shape != expected_shape: @@ -212,6 +250,9 @@ def _section_branches( if section.backward_modes is None: backward_modes = [(e_field.copy(), -h_field.copy()) for e_field, h_field in forward_modes] backward_wavenumbers = -wavenumbers + if section.backward_wavenumber_matrix is not None: + raise ValueError('backward_wavenumber_matrix requires explicit backward modes and wavenumbers') + backward_wavenumber_matrix = -wavenumber_matrix else: backward_shape, _backward_h_shape = _validate_port_modes( f'sections[{index}].backward_modes', @@ -225,6 +266,11 @@ def _section_branches( section.backward_wavenumbers, ) backward_modes = _as_mode_arrays(section.backward_modes) + backward_wavenumber_matrix = _as_wavenumber_matrix( + f'sections[{index}].backward_wavenumber_matrix', + section.backward_wavenumber_matrix, + backward_wavenumbers, + ) if len(backward_wavenumbers) != len(wavenumbers): raise ValueError('forward and backward mode counts must match') @@ -259,7 +305,10 @@ def _section_branches( if len(dual_modes) != len(forward_modes) or len(dual_backward_modes) != len(backward_modes): raise ValueError('dual mode counts must match the corresponding modal basis counts') - return z_coord, forward_modes, wavenumbers, backward_modes, backward_wavenumbers, dual_modes, dual_backward_modes, shape + return ( + z_coord, forward_modes, wavenumbers, backward_modes, backward_wavenumbers, + dual_modes, dual_backward_modes, wavenumber_matrix, backward_wavenumber_matrix, shape, + ) def _validate_taper_sections( @@ -272,6 +321,7 @@ def _validate_taper_sections( list[NDArray[numpy.complex128]], int, ]: + """Validate, phase-align, and assemble all taper branch data.""" if len(sections) < 2: raise ValueError('at least two taper sections are required') @@ -279,12 +329,11 @@ def _validate_taper_sections( branch_modes: list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]] = [] branch_dual_modes: list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]] = [] branch_wavenumbers: list[NDArray[numpy.complex128]] = [] - explicit_duals: list[bool] = [] expected_count: int | None = None expected_shape: tuple[int, ...] | None = None for index, section in enumerate(sections): - z_coord, forward_modes, forward_wavenumbers, backward_modes, backward_wavenumbers, dual_modes, dual_backward_modes, shape = _section_branches( + z_coord, forward_modes, forward_wavenumbers, backward_modes, backward_wavenumbers, dual_modes, dual_backward_modes, forward_matrix, backward_matrix, shape = _section_branches( section, index, expected_count, @@ -296,17 +345,22 @@ def _validate_taper_sections( z_coords.append(z_coord) branch_modes.append([*forward_modes, *backward_modes]) branch_dual_modes.append([*dual_modes, *dual_backward_modes]) - branch_wavenumbers.append(numpy.concatenate((forward_wavenumbers, backward_wavenumbers))) - explicit_duals.append(section.dual_modes is not None or section.dual_backward_modes is not None) + branch_wavenumbers.append(linalg.block_diag(forward_matrix, backward_matrix)) z_array = numpy.asarray(z_coords, dtype=float) if not (numpy.diff(z_array) > 0).all(): raise ValueError('taper section z coordinates must be strictly increasing') for index in range(1, len(branch_modes)): - branch_modes[index] = _phase_align_modes(branch_modes[index - 1], branch_modes[index], dxes, branch_dual_modes[index - 1]) - if not explicit_duals[index]: - branch_dual_modes[index] = branch_modes[index] + branch_modes[index], transform = _phase_align_modes( + branch_modes[index - 1], branch_modes[index], dxes, branch_dual_modes[index - 1], + ) + inverse = numpy.linalg.inv(transform) + branch_dual_modes[index] = [ + (e_field * inverse[row, row], h_field * inverse[row, row]) + for row, (e_field, h_field) in enumerate(branch_dual_modes[index]) + ] + branch_wavenumbers[index] = inverse @ branch_wavenumbers[index] @ transform assert expected_count is not None return z_array, branch_modes, branch_dual_modes, branch_wavenumbers, expected_count @@ -321,6 +375,7 @@ def _taper_interval_generator( dz: float, dxes: dx_lists2_t, ) -> NDArray[numpy.complex128]: + """Build the local coupled-mode generator for one sampled interval.""" mode_count = len(left_modes) gram = numpy.zeros((mode_count, mode_count), dtype=complex) derivative_overlap = numpy.zeros((mode_count, mode_count), dtype=complex) @@ -336,7 +391,7 @@ def _taper_interval_generator( derivative_overlap[row, col] = _lorentz_derivative_overlap(left_row_mode, derivative, dxes) coupling = numpy.linalg.pinv(gram) @ derivative_overlap - propagation = numpy.diag(-1j * 0.5 * (left_wavenumbers + right_wavenumbers)) + propagation = -1j * 0.5 * (left_wavenumbers + right_wavenumbers) return propagation - coupling @@ -344,6 +399,7 @@ def _abcd_to_s( abcd: NDArray[numpy.complex128], n_modes: int, ) -> NDArray[numpy.complex128]: + """Convert a block transfer matrix to the corresponding scattering matrix.""" A = abcd[:n_modes, :n_modes] B = abcd[:n_modes, n_modes:] C = abcd[n_modes:, :n_modes] @@ -586,3 +642,123 @@ def get_taper_s( ss = 0.5 * (ss + ss.T) return ss + + +def mode_sections_from_bases( + z_coords: Sequence[float], + bases: Sequence[ComplexModeBasis], + *, + mode_indices: Sequence[int] | None = None, + backward_bases: Sequence[ComplexModeBasis] | None = None, + coupling_rtol: float = 1e-10, + ) -> tuple[ModeSection, ...]: + """Convert tracked complex-mode bases into EME taper sections. + + The caller is responsible for tracking the bases first and for explicitly + power-normalizing any port modes whose scattering amplitudes will be + interpreted as power waves. When ``backward_bases`` is omitted, the normal + reciprocal ``(E, -H, -beta)`` branches are synthesized by ``ModeSection``. + A selected mode subset must be invariant under each basis's wavenumber + matrix; ``coupling_rtol`` controls the relative numerical coupling threshold. + + Args: + z_coords: Strictly increasing longitudinal coordinates for the bases. + bases: Equal-sized, already tracked forward modal bases. + mode_indices: Optional ordered subset of modal basis indices. Omit to + retain every mode. + backward_bases: Optional explicitly solved backward bases. When + omitted, reciprocal backward branches are synthesized later. + coupling_rtol: Relative threshold below which propagation-matrix + coupling across the selected subset is treated as numerical zero. + + Returns: + Immutable sequence of EME mode sections carrying primal, dual, and + full propagation-matrix data. + + Raises: + ValueError: If sequence lengths, mode indices, tolerances, or selected + propagation blocks are incompatible. + """ + if len(z_coords) != len(bases): + raise ValueError('z_coords and bases must have the same length') + if len(bases) < 2: + raise ValueError('at least two bases are required') + if backward_bases is not None and len(backward_bases) != len(bases): + raise ValueError('backward_bases must have the same length as bases') + if not numpy.isfinite(coupling_rtol) or coupling_rtol <= 0: + raise ValueError('coupling_rtol must be finite and positive') + + if mode_indices is None: + indices = numpy.arange(len(bases[0]), dtype=int) + else: + raw_indices = numpy.asarray(mode_indices) + if raw_indices.ndim != 1 or raw_indices.size == 0: + raise ValueError('mode_indices must be a nonempty one-dimensional sequence') + if any(isinstance(index, bool | numpy.bool_) for index in mode_indices): + raise ValueError('mode_indices entries must be integers') + try: + indices = numpy.asarray([operator.index(index) for index in mode_indices], dtype=int) + except TypeError as error: + raise ValueError('mode_indices entries must be integers') from error + if indices.ndim != 1 or indices.size == 0: + raise ValueError('mode_indices must be a nonempty one-dimensional sequence') + if numpy.any(indices < 0): + raise ValueError('mode_indices must be nonnegative') + if numpy.unique(indices).size != indices.size: + raise ValueError('mode_indices must not contain duplicates') + + def validate_invariant_subset(matrix: NDArray[numpy.complex128], size: int, section_name: str) -> None: + """Reject a subset coupled to excluded coordinates by its beta matrix.""" + excluded = numpy.setdiff1d(numpy.arange(size), indices, assume_unique=False) + if excluded.size == 0: + return + scale = max(1.0, float(numpy.linalg.norm(matrix))) + coupling = max( + float(numpy.linalg.norm(matrix[numpy.ix_(indices, excluded)])), + float(numpy.linalg.norm(matrix[numpy.ix_(excluded, indices)])), + ) + if coupling > coupling_rtol * scale: + raise ValueError( + f'mode_indices cut through a coupled wavenumber block in {section_name}; ' + 'select the complete coupled block', + ) + + sections: list[ModeSection] = [] + for section_index, (z_coord, basis) in enumerate(zip(z_coords, bases, strict=True)): + if numpy.any(indices >= len(basis)): + raise ValueError(f'mode_indices are out of range for bases[{section_index}]') + validate_invariant_subset(basis.wavenumber_matrix, len(basis), f'bases[{section_index}]') + modes = [basis.modes[index] for index in indices] + dual_modes = [basis.dual_modes[index] for index in indices] + wavenumbers = basis.wavenumbers[indices] + wavenumber_matrix = basis.wavenumber_matrix[numpy.ix_(indices, indices)] + if backward_bases is None: + sections.append(ModeSection( + z=float(z_coord), + modes=modes, + wavenumbers=wavenumbers, + dual_modes=dual_modes, + wavenumber_matrix=wavenumber_matrix, + )) + continue + + backward_basis = backward_bases[section_index] + if numpy.any(indices >= len(backward_basis)): + raise ValueError(f'mode_indices are out of range for backward_bases[{section_index}]') + validate_invariant_subset( + backward_basis.wavenumber_matrix, + len(backward_basis), + f'backward_bases[{section_index}]', + ) + sections.append(ModeSection( + z=float(z_coord), + modes=modes, + wavenumbers=wavenumbers, + backward_modes=[backward_basis.modes[index] for index in indices], + backward_wavenumbers=backward_basis.wavenumbers[indices], + dual_modes=dual_modes, + dual_backward_modes=[backward_basis.dual_modes[index] for index in indices], + wavenumber_matrix=wavenumber_matrix, + backward_wavenumber_matrix=backward_basis.wavenumber_matrix[numpy.ix_(indices, indices)], + )) + return tuple(sections) diff --git a/meanas/fdfd/waveguide_2d_complex.py b/meanas/fdfd/waveguide_2d_complex.py new file mode 100644 index 0000000..eadf316 --- /dev/null +++ b/meanas/fdfd/waveguide_2d_complex.py @@ -0,0 +1,881 @@ +"""Targeted complex waveguide modes on open, complex-stretched cross-sections.""" +from collections.abc import Callable, Sequence +import dataclasses +from typing import TYPE_CHECKING, Any + +import numpy +from scipy import sparse +from scipy import linalg +from scipy.optimize import linear_sum_assignment +from scipy.sparse import linalg as splinalg + +from ..eigensolvers import TargetedEigensolution, scipy_eigs_near +from ..fdmath.operators import deriv_back, deriv_forward +from . import waveguide_2d + +if TYPE_CHECKING: + from numpy.typing import NDArray + + from ..fdmath import dx_lists2_t, vcfdfield2, vfdslice + + +type targeted_eigensolver_t = Callable[..., TargetedEigensolution] + + +@dataclasses.dataclass(frozen=True) +class ComplexModeDiagnostics: + """Per-mode numerical and physical diagnostics. + + ``biorthogonality_condition`` is the condition number of the right/left + Gram matrix after independently equilibrating its vector columns. + + Attributes: + right_residuals: Relative residual of each primal mode or basis column. + left_residuals: Relative residual of each transpose-left mode or column. + powers: Conjugated Poynting power of each primal modal field. + pml_participation: Fraction of transverse state weight inside the PML. + target_indices: Index of the nearest requested target for each spectral + value. + target_distances: Absolute distance from each spectral value to that + target. + biorthogonality_error: Maximum error in the Lorentz dual Gram matrix. + biorthogonality_condition: Equilibrated right/left Gram condition + number. + """ + + right_residuals: NDArray[numpy.float64] + left_residuals: NDArray[numpy.float64] + powers: NDArray[numpy.complex128] + pml_participation: NDArray[numpy.float64] + target_indices: NDArray[numpy.int_] + target_distances: NDArray[numpy.float64] + biorthogonality_error: float + biorthogonality_condition: float + + +@dataclasses.dataclass(frozen=True) +class ComplexModeBasis: + """Primal and Lorentz-dual modes returned by :func:`solve_modes_near`. + + ``wavenumbers`` records the eigenspectrum. ``wavenumber_matrix`` is the + actual beta operator in the current field basis; it becomes non-diagonal + when tracking rotates a nearly degenerate subspace. + + Attributes: + wavenumbers: Eigenvalue spectrum associated with the basis. + wavenumber_matrix: Representation of the beta operator in the current + possibly mixed basis. + electric_fields: Primal vectorized three-component electric fields. + magnetic_fields: Primal vectorized three-component magnetic fields. + dual_electric_fields: Lorentz-dual electric projection fields. + dual_magnetic_fields: Lorentz-dual magnetic projection fields. + diagnostics: Numerical and physical diagnostics for the basis. + solve_dxes: Possibly complex-stretched grid used by the eigenproblem. + overlap_dxes: Real-space grid used for overlaps and normalization. + omega: Angular frequency used to construct the operator. + epsilon: Vectorized three-component relative permittivity. + mu: Optional vectorized relative permeability. + physical_mask: Transverse cells treated as outside the PML. + """ + + wavenumbers: NDArray[numpy.complex128] + wavenumber_matrix: NDArray[numpy.complex128] + electric_fields: NDArray[numpy.complex128] + magnetic_fields: NDArray[numpy.complex128] + dual_electric_fields: NDArray[numpy.complex128] + dual_magnetic_fields: NDArray[numpy.complex128] + diagnostics: ComplexModeDiagnostics + solve_dxes: dx_lists2_t + overlap_dxes: dx_lists2_t + omega: complex + epsilon: NDArray[numpy.complex128] + mu: NDArray[numpy.complex128] | None + physical_mask: NDArray[numpy.bool_] + + def __len__(self) -> int: + """Return the number of modal basis vectors.""" + return self.wavenumbers.size + + @property + def modes(self) -> tuple[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]], ...]: + """Return primal electric/magnetic field pairs in basis order.""" + return tuple(zip(self.electric_fields, self.magnetic_fields, strict=True)) + + @property + def dual_modes(self) -> tuple[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]], ...]: + """Return Lorentz-dual electric/magnetic field pairs in basis order.""" + return tuple(zip(self.dual_electric_fields, self.dual_magnetic_fields, strict=True)) + + +def operator_beta( + omega: complex, + dxes: dx_lists2_t, + epsilon: vfdslice, + mu: vfdslice | None = None, + ) -> sparse.sparray: + r"""Build the first-order transverse Maxwell operator ``M psi = beta psi``. + + The state ordering is ``[E_x, E_y, H_x, H_y]``. Longitudinal fields are + eliminated, but both transverse field families are retained so the sign of + ``beta`` remains explicit. Algebraically, ``M @ M`` has the existing + electric and magnetic ``beta**2`` waveguide operators on its diagonal. + + Args: + omega: Nonzero angular frequency. + dxes: Two-dimensional Yee-grid spacings, which may include complex PML + stretching. + epsilon: Vectorized three-component relative permittivity. + mu: Optional vectorized three-component relative permeability. Unit + permeability is used when omitted. + + Returns: + Sparse first-order signed-beta operator on transverse E/H states. + + Raises: + ValueError: If the material vectors, grid shape, or frequency are + incompatible. + """ + epsilon = numpy.asarray(epsilon) + if mu is None: + mu = numpy.ones_like(epsilon) + mu = numpy.asarray(mu) + if epsilon.ndim != 1 or mu.ndim != 1 or epsilon.shape != mu.shape or epsilon.size % 3: + raise ValueError('epsilon and mu must be equal-length vectorized three-component fields') + if omega == 0: + raise ValueError('omega must be nonzero') + + n_points = epsilon.size // 3 + if n_points != dxes[0][0].size * dxes[0][1].size: + raise ValueError('material and grid shapes do not match') + + dfx, dfy = deriv_forward(dxes[0]) + dbx, dby = deriv_back(dxes[1]) + zero = sparse.csr_array((n_points, n_points)) + identity = sparse.eye_array(n_points) + rotation = sparse.block_array([[zero, -identity], [identity, zero]], format='csr') + + eps_x, eps_y, eps_z = numpy.split(epsilon, 3) + mu_x, mu_y, mu_z = numpy.split(mu, 3) + eps_xy = sparse.diags_array(numpy.concatenate((eps_x, eps_y))) + mu_xy = sparse.diags_array(numpy.concatenate((mu_x, mu_y))) + eps_z_inv = sparse.diags_array(1 / eps_z) + mu_z_inv = sparse.diags_array(1 / mu_z) + + grad_forward = sparse.vstack((dfx, dfy), format='csr') + grad_backward = sparse.vstack((dbx, dby), format='csr') + rotated_div_backward = sparse.hstack((-dby, dbx), format='csr') + rotated_div_forward = sparse.hstack((-dfy, dfx), format='csr') + + beta_e_from_h = ( + -omega * rotation @ mu_xy + + grad_forward @ eps_z_inv @ rotated_div_backward / omega + ) + beta_h_from_e = ( + omega * rotation @ eps_xy + - grad_backward @ mu_z_inv @ rotated_div_forward / omega + ) + empty = sparse.csr_array((2 * n_points, 2 * n_points)) + return sparse.block_array( + [[empty, beta_e_from_h], [beta_h_from_e, empty]], + format = 'csr', + ) + + +def lorentz_metric(dxes: dx_lists2_t) -> sparse.sparray: + """Return the state metric matching the EME symmetric Lorentz overlap. + + Args: + dxes: Two-dimensional Yee-grid spacings used for physical overlaps. + + Returns: + Sparse bilinear metric on states ordered as ``[Ex, Ey, Hx, Hy]``. + """ + nx = dxes[0][0].size + ny = dxes[0][1].size + real_dxes = [[numpy.real(values) for values in family] for family in dxes] + weight_ex_hy = numpy.outer(real_dxes[1][0], real_dxes[0][1]).reshape(nx * ny) + weight_ey_hx = numpy.outer(real_dxes[0][0], real_dxes[1][1]).reshape(nx * ny) + diagonal_a = sparse.diags_array(weight_ex_hy) + diagonal_b = sparse.diags_array(weight_ey_hx) + empty_n = sparse.csr_array((nx * ny, nx * ny)) + # inner_product(E, H) = E.T @ q @ H; it already includes a factor 1/2. + q_matrix = 0.5 * sparse.block_array( + [[empty_n, diagonal_a], [-diagonal_b, empty_n]], + format = 'csr', + ) + empty_2n = sparse.csr_array(q_matrix.shape) + # _lorentz_overlap adds the two reaction products and multiplies by 1/2. + return 0.5 * sparse.block_array( + [[empty_2n, q_matrix], [q_matrix.T, empty_2n]], + format = 'csr', + ) + + +def _copy_dxes(dxes: dx_lists2_t) -> list[list[NDArray]]: + """Deep-copy the two Yee-grid spacing families.""" + return [[numpy.asarray(values).copy() for values in family] for family in dxes] + + +def _infer_physical_mask(dxes: dx_lists2_t) -> NDArray[numpy.bool_]: + """Infer the non-PML aperture from imaginary grid stretching.""" + pml_x = numpy.zeros(dxes[0][0].size, dtype=bool) + pml_y = numpy.zeros(dxes[0][1].size, dtype=bool) + for family in dxes: + pml_x |= numpy.abs(numpy.imag(family[0])) > 0 + pml_y |= numpy.abs(numpy.imag(family[1])) > 0 + return ~(pml_x[:, None] | pml_y[None, :]) + + +def _state_weights(dxes: dx_lists2_t, physical_mask: NDArray[numpy.bool_]) -> NDArray[numpy.float64]: + """Build cell-area weights for all four transverse state components.""" + dx = 0.5 * (numpy.real(dxes[0][0]) + numpy.real(dxes[1][0])) + dy = 0.5 * (numpy.real(dxes[0][1]) + numpy.real(dxes[1][1])) + weights = numpy.outer(dx, dy) * physical_mask + return numpy.tile(weights.ravel(), 4) + + +def _normalize_columns( + matrix: NDArray[numpy.complex128], + ) -> tuple[NDArray[numpy.complex128], NDArray[numpy.float64]]: + """Normalize matrix columns and return their original Euclidean norms.""" + norms = numpy.linalg.norm(matrix, axis=0) + normalized = numpy.divide(matrix, norms, out=numpy.zeros_like(matrix), where=norms != 0) + return normalized, numpy.asarray(norms, dtype=float) + + +def _select_unique_modes( + solutions: Sequence[TargetedEigensolution], + targets: NDArray[numpy.complex128], + modes_per_target: int, + *, + eigenvalue_rtol: float, + subspace_tol: float, + state_weights: NDArray[numpy.float64], + ) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128], NDArray[numpy.complex128], NDArray[numpy.int_]]: + """Merge targeted candidates into independent, cross-pairable eigenspaces.""" + sqrt_weights = numpy.sqrt(state_weights) + + def rank_and_pivots(vectors: list[NDArray[numpy.complex128]]) -> tuple[int, NDArray[numpy.int_]]: + """Estimate weighted subspace rank and QR pivot order.""" + if not vectors: + return 0, numpy.empty(0, dtype=int) + matrix = sqrt_weights[:, None] * numpy.column_stack(vectors) + matrix, _norms = _normalize_columns(matrix) + _q, triangular, pivots = linalg.qr(matrix, mode='economic', pivoting=True) + diagonal = numpy.abs(numpy.diag(triangular)) + if diagonal.size == 0 or diagonal[0] == 0: + return 0, numpy.asarray(pivots, dtype=int) + rank = int(numpy.count_nonzero(diagonal > subspace_tol * diagonal[0])) + return rank, numpy.asarray(pivots, dtype=int) + + def normalized_cross_gram( + right_vectors: list[NDArray[numpy.complex128]], + left_vectors: list[NDArray[numpy.complex128]], + ) -> NDArray[numpy.complex128]: + """Form the bilinear cross Gram matrix after column equilibration.""" + right = numpy.column_stack(right_vectors) + left = numpy.column_stack(left_vectors) + right, _right_norms = _normalize_columns(right) + left, _left_norms = _normalize_columns(left) + return left.T @ right + + def matrix_rank(matrix: NDArray[numpy.complex128]) -> int: + """Estimate matrix rank using the configured subspace tolerance.""" + singular_values = linalg.svdvals(matrix) + if singular_values.size == 0: + return 0 + cutoff = subspace_tol * max(1.0, float(singular_values[0])) + return int(numpy.count_nonzero(singular_values > cutoff)) + + initial_records: list[tuple[complex, NDArray[numpy.complex128], NDArray[numpy.complex128]]] = [] + remaining_records: list[tuple[complex, NDArray[numpy.complex128], NDArray[numpy.complex128]]] = [] + for target_index, (target, solution) in enumerate(zip(targets, solutions, strict=True)): + local_right: list[NDArray[numpy.complex128]] = [] + local_left: list[NDArray[numpy.complex128]] = [] + order = numpy.argsort(numpy.abs(solution.eigenvalues - target)) + stop_position: int | None = None + for position, candidate_index in enumerate(order): + value = complex(solution.eigenvalues[candidate_index]) + right_vector = solution.right_eigenvectors[:, candidate_index] + left_vector = solution.left_eigenvectors[:, candidate_index] + initial_records.append((value, right_vector, left_vector)) + local_right.append(right_vector) + local_left.append(left_vector) + right_rank, _right_pivots = rank_and_pivots(local_right) + left_rank, _left_pivots = rank_and_pivots(local_left) + if right_rank >= modes_per_target and left_rank >= modes_per_target: + stop_position = position + break + if stop_position is None: + raise RuntimeError(f'target {target_index} did not produce enough independent right and left modes') + for candidate_index in order[stop_position + 1:]: + remaining_records.append(( + complex(solution.eigenvalues[candidate_index]), + solution.right_eigenvectors[:, candidate_index], + solution.left_eigenvectors[:, candidate_index], + )) + + clusters: list[list[int]] = [] + for candidate_index, (value, _right, _left) in enumerate(initial_records): + for cluster in clusters: + anchor = initial_records[cluster[0]][0] + if abs(value - anchor) <= eigenvalue_rtol * max(1.0, abs(value), abs(anchor)): + cluster.append(candidate_index) + break + else: + clusters.append([candidate_index]) + + chosen_values: list[complex] = [] + chosen_right: list[NDArray[numpy.complex128]] = [] + chosen_left: list[NDArray[numpy.complex128]] = [] + for cluster in clusters: + anchor = initial_records[cluster[0]][0] + right_values = [initial_records[index][0] for index in cluster] + left_values = [initial_records[index][0] for index in cluster] + right_vectors = [initial_records[index][1] for index in cluster] + left_vectors = [initial_records[index][2] for index in cluster] + right_rank, _right_pivots = rank_and_pivots(right_vectors) + left_rank, _left_pivots = rank_and_pivots(left_vectors) + required_rank = max(right_rank, left_rank) + cross_gram = normalized_cross_gram(right_vectors, left_vectors) + cross_rank = matrix_rank(cross_gram) + for value, right_vector, left_vector in remaining_records: + if cross_rank >= required_rank: + break + if abs(value - anchor) > eigenvalue_rtol * max(1.0, abs(value), abs(anchor)): + continue + right_vectors.append(right_vector) + left_vectors.append(left_vector) + right_values.append(value) + left_values.append(value) + cross_gram = normalized_cross_gram(right_vectors, left_vectors) + cross_rank = matrix_rank(cross_gram) + if cross_rank < required_rank: + raise RuntimeError( + f'eigenvalue cluster near {anchor!r} has incompatible right/transpose-left ' + f'subspaces (required cross rank {required_rank}, available {cross_rank})', + ) + + _q, _triangular, right_pivots = linalg.qr(cross_gram, mode='economic', pivoting=True) + right_selection = numpy.sort(numpy.asarray(right_pivots[:required_rank], dtype=int)) + _q, _triangular, left_pivots = linalg.qr( + cross_gram[:, right_selection].T, + mode='economic', + pivoting=True, + ) + left_selection = numpy.asarray(left_pivots[:required_rank], dtype=int) + selected_cross_gram = cross_gram[numpy.ix_(left_selection, right_selection)] + selected_cross_rank = matrix_rank(selected_cross_gram) + if selected_cross_rank < required_rank: + raise RuntimeError( + f'eigenvalue cluster near {anchor!r} has incompatible right/transpose-left ' + f'subspaces (required cross rank {required_rank}, selected {selected_cross_rank})', + ) + selected_right_values = numpy.asarray(right_values)[right_selection] + selected_left_values = numpy.asarray(left_values)[left_selection] + rows, columns = linear_sum_assignment(numpy.abs(selected_right_values[:, None] - selected_left_values[None, :])) + paired_left = numpy.asarray(left_selection)[columns[numpy.argsort(rows)]] + chosen_values.extend(numpy.asarray(right_values)[right_selection]) + chosen_right.extend(right_vectors[int(index)] for index in right_selection) + chosen_left.extend(left_vectors[int(index)] for index in paired_left) + + values = numpy.asarray(chosen_values, dtype=complex) + target_indices = numpy.argmin(numpy.abs(values[:, None] - targets[None, :]), axis=1) + return ( + values, + numpy.column_stack(chosen_right), + numpy.column_stack(chosen_left), + numpy.asarray(target_indices, dtype=int), + ) + + +def _relative_residuals( + operator: sparse.sparray, + values: NDArray[numpy.complex128], + vectors: NDArray[numpy.complex128], + ) -> NDArray[numpy.float64]: + """Compute scale-aware eigenpair residuals for individual vectors.""" + residuals = numpy.empty(values.size, dtype=float) + for index, value in enumerate(values): + vector = vectors[:, index] + residuals[index] = numpy.linalg.norm(operator @ vector - value * vector) / ( + max(1.0, abs(value)) * numpy.linalg.norm(vector) + ) + return residuals + + +def solve_modes_near( + omega: complex, + solve_dxes: dx_lists2_t, + epsilon: vfdslice, + *, + target_wavenumbers: Sequence[complex], + modes_per_target: int, + mu: vfdslice | None = None, + overlap_dxes: dx_lists2_t | None = None, + physical_mask: NDArray[numpy.bool_] | None = None, + candidate_margin: int = 4, + eigensolver: targeted_eigensolver_t = scipy_eigs_near, + eigensolver_opts: dict[str, Any] | None = None, + residual_tol: float = 1e-8, + duplicate_eigenvalue_rtol: float = 1e-8, + duplicate_subspace_tol: float = 1e-7, + ) -> ComplexModeBasis: + """Solve for complex signed waveguide modes near one or more beta targets. + + Each target requests an independent candidate solve. Overlapping results + are merged by eigenvalue cluster and weighted subspace rank, after which + transpose-left eigenvectors are converted into a Lorentz-biorthogonal dual + basis. Returned fields use the full three-component vectorization, while + the underlying eigenproblem is ordered as ``[Ex, Ey, Hx, Hy]``. + + Args: + omega: Nonzero angular frequency. + solve_dxes: Yee-grid spacings used to build the possibly PML-stretched + eigenproblem. + epsilon: Vectorized three-component relative permittivity. + target_wavenumbers: Finite signed-beta targets for shift-invert solves. + modes_per_target: Independent candidates requested around each target. + mu: Optional vectorized three-component relative permeability. + overlap_dxes: Grid used for physical overlaps and normalization. By + default, the real parts of ``solve_dxes`` are used. + physical_mask: Optional transverse mask selecting the non-PML aperture. + By default, cells with imaginary grid stretching are excluded. + candidate_margin: Extra eigenpairs requested from each targeted solve + to support deduplication and degenerate-subspace repair. + eigensolver: Backend returning paired right and transpose-left vectors. + eigensolver_opts: Keyword arguments forwarded to ``eigensolver``. + residual_tol: Maximum accepted right, left, and dual-basis residual. + duplicate_eigenvalue_rtol: Relative tolerance for eigenvalue clustering. + duplicate_subspace_tol: Relative numerical-rank tolerance used while + merging repeated or degenerate candidates. + + Returns: + Complex modal basis with primal fields, Lorentz duals, propagation data, + copied solve metadata, and numerical diagnostics. + + Raises: + ValueError: If inputs or backend result shapes contain invalid values. + RuntimeError: If too few independent modes are found, right/left spaces + cannot be paired, or requested residual tolerances are not met. + """ + if modes_per_target < 1: + raise ValueError('modes_per_target must be positive') + if candidate_margin < 0: + raise ValueError('candidate_margin must be nonnegative') + if not numpy.isfinite(residual_tol) or residual_tol <= 0: + raise ValueError('residual_tol must be finite and positive') + if not numpy.isfinite(duplicate_eigenvalue_rtol) or duplicate_eigenvalue_rtol <= 0: + raise ValueError('duplicate_eigenvalue_rtol must be finite and positive') + if not numpy.isfinite(duplicate_subspace_tol) or duplicate_subspace_tol <= 0: + raise ValueError('duplicate_subspace_tol must be finite and positive') + targets = numpy.asarray(target_wavenumbers, dtype=complex) + if targets.ndim != 1 or targets.size == 0 or not numpy.isfinite(targets).all(): + raise ValueError('target_wavenumbers must be a finite nonempty one-dimensional sequence') + if overlap_dxes is None: + overlap_dxes = [[numpy.real(values) for values in family] for family in solve_dxes] + if physical_mask is None: + physical_mask = _infer_physical_mask(solve_dxes) + physical_mask = numpy.asarray(physical_mask, dtype=bool) + grid_shape = (solve_dxes[0][0].size, solve_dxes[0][1].size) + if physical_mask.shape != grid_shape: + raise ValueError('physical_mask must match the transverse grid shape') + + operator = operator_beta(omega, solve_dxes, epsilon, mu) + options = {} if eigensolver_opts is None else dict(eigensolver_opts) + solutions = [ + eigensolver(operator, complex(target), modes_per_target + candidate_margin, **options) + for target in targets + ] + for target_index, solution in enumerate(solutions): + count = numpy.asarray(solution.eigenvalues).size + if numpy.asarray(solution.eigenvalues).shape != (count,): + raise ValueError(f'eigensolver result {target_index} eigenvalues must be one-dimensional') + expected_shape = (operator.shape[0], count) + if numpy.asarray(solution.right_eigenvectors).shape != expected_shape or numpy.asarray(solution.left_eigenvectors).shape != expected_shape: + raise ValueError(f'eigensolver result {target_index} eigenvector shapes are invalid') + if count < modes_per_target: + raise RuntimeError(f'eigensolver result {target_index} returned too few candidates') + if not ( + numpy.isfinite(solution.eigenvalues).all() + and numpy.isfinite(solution.right_eigenvectors).all() + and numpy.isfinite(solution.left_eigenvectors).all() + ): + raise ValueError(f'eigensolver result {target_index} contains non-finite values') + merge_weights = _state_weights(overlap_dxes, numpy.ones(grid_shape, dtype=bool)) + values, right, left, target_indices = _select_unique_modes( + solutions, + targets, + modes_per_target, + eigenvalue_rtol = duplicate_eigenvalue_rtol, + subspace_tol = duplicate_subspace_tol, + state_weights = merge_weights, + ) + + right, right_norms = _normalize_columns(right) + left, left_norms = _normalize_columns(left) + if numpy.any(right_norms == 0) or numpy.any(left_norms == 0): + raise RuntimeError('selected right/transpose-left eigenbasis contains a zero vector') + gram = left.T @ right + condition = float(numpy.linalg.cond(gram)) + if not numpy.isfinite(condition) or condition > 1e12: + raise RuntimeError(f'left/right eigenbasis is singular or ill-conditioned ({condition=:.3e})') + left = linalg.solve(gram, left.T).T + + weights = _state_weights(overlap_dxes, physical_mask) + for index in range(values.size): + vector = right[:, index] + norm_value = numpy.sqrt(numpy.sum(weights * numpy.abs(vector) ** 2)) + if norm_value == 0: + raise RuntimeError(f'mode {index} has no support in the physical aperture') + weighted_vector = numpy.where(weights > 0, vector, 0) + phase_sample = weighted_vector[numpy.argmax(numpy.abs(weighted_vector))] + scale = numpy.exp(-1j * numpy.angle(phase_sample)) / norm_value + right[:, index] *= scale + left[:, index] /= scale + + right_residuals = _relative_residuals(operator, values, right) + left_residuals = _relative_residuals(operator.T, values, left) + if numpy.max(right_residuals) > residual_tol or numpy.max(left_residuals) > residual_tol: + raise RuntimeError('targeted eigenmodes did not meet the requested residual tolerance') + + metric = lorentz_metric(overlap_dxes).tocsc() + dual_states = numpy.asarray(splinalg.spsolve(metric.T, left)) + if dual_states.ndim == 1: + dual_states = dual_states[:, None] + biorthogonal_gram = dual_states.T @ metric @ right + biorthogonality_error = float(numpy.max(numpy.abs(biorthogonal_gram - numpy.eye(values.size)))) + if biorthogonality_error > residual_tol: + raise RuntimeError('failed to construct a biorthogonal Lorentz dual basis') + + n_points = grid_shape[0] * grid_shape[1] + electric_fields = numpy.empty((values.size, 3 * n_points), dtype=complex) + magnetic_fields = numpy.empty_like(electric_fields) + dual_electric_fields = numpy.zeros_like(electric_fields) + dual_magnetic_fields = numpy.zeros_like(electric_fields) + epsilon_array = numpy.asarray(epsilon) + mu_array = None if mu is None else numpy.asarray(mu) + for index, beta in enumerate(values): + e_xy = right[:2 * n_points, index] + h_xy = right[2 * n_points:, index] + electric_fields[index] = waveguide_2d.exy2e(beta, solve_dxes, epsilon_array) @ e_xy + magnetic_fields[index] = waveguide_2d.hxy2h(beta, solve_dxes, mu_array) @ h_xy + dual_electric_fields[index, :2 * n_points] = dual_states[:2 * n_points, index] + dual_magnetic_fields[index, :2 * n_points] = dual_states[2 * n_points:, index] + + powers = numpy.asarray([ + waveguide_2d.inner_product(e_field, h_field, overlap_dxes, conj_h=True) + for e_field, h_field in zip(electric_fields, magnetic_fields, strict=True) + ]) + full_weights = numpy.tile(numpy.outer( + 0.5 * (numpy.real(overlap_dxes[0][0]) + numpy.real(overlap_dxes[1][0])), + 0.5 * (numpy.real(overlap_dxes[0][1]) + numpy.real(overlap_dxes[1][1])), + ).ravel(), 4) + pml_weights = full_weights * numpy.tile(~physical_mask.ravel(), 4) + pml_participation = numpy.sum(pml_weights[:, None] * numpy.abs(right) ** 2, axis=0) / numpy.sum( + full_weights[:, None] * numpy.abs(right) ** 2, + axis = 0, + ) + target_distances = numpy.abs(values - targets[target_indices]) + + diagnostics = ComplexModeDiagnostics( + right_residuals = right_residuals, + left_residuals = left_residuals, + powers = powers, + pml_participation = numpy.asarray(pml_participation, dtype=float), + target_indices = target_indices, + target_distances = numpy.asarray(target_distances, dtype=float), + biorthogonality_error = biorthogonality_error, + biorthogonality_condition = condition, + ) + return ComplexModeBasis( + wavenumbers = values, + wavenumber_matrix = numpy.diag(values), + electric_fields = electric_fields, + magnetic_fields = magnetic_fields, + dual_electric_fields = dual_electric_fields, + dual_magnetic_fields = dual_magnetic_fields, + diagnostics = diagnostics, + solve_dxes = _copy_dxes(solve_dxes), + overlap_dxes = _copy_dxes(overlap_dxes), + omega = complex(omega), + epsilon = numpy.asarray(epsilon, dtype=complex).copy(), + mu = None if mu is None else numpy.asarray(mu, dtype=complex).copy(), + physical_mask = physical_mask.copy(), + ) + + +def _block_relative_residuals( + operator: sparse.sparray, + matrix: NDArray[numpy.complex128], + vectors: NDArray[numpy.complex128], + ) -> NDArray[numpy.float64]: + """Compute column residuals for an invariant subspace representation.""" + image = operator @ vectors + expected = vectors @ matrix + numerator = numpy.linalg.norm(image - expected, axis=0) + denominator = numpy.maximum(numpy.linalg.norm(image, axis=0), numpy.linalg.norm(expected, axis=0)) + return numpy.asarray(numerator / numpy.maximum(denominator, numpy.finfo(float).tiny), dtype=float) + + +def _refresh_diagnostics(basis: ComplexModeBasis) -> ComplexModeBasis: + """Recompute basis-dependent residual, power, PML, and Gram diagnostics.""" + n_points = basis.electric_fields.shape[1] // 3 + right = numpy.vstack((basis.electric_fields[:, :2 * n_points].T, basis.magnetic_fields[:, :2 * n_points].T)) + dual = numpy.vstack((basis.dual_electric_fields[:, :2 * n_points].T, basis.dual_magnetic_fields[:, :2 * n_points].T)) + operator = operator_beta(basis.omega, basis.solve_dxes, basis.epsilon, basis.mu) + metric = lorentz_metric(basis.overlap_dxes) + left = metric.T @ dual + gram = dual.T @ metric @ right + normalized_right, right_norms = _normalize_columns(right) + normalized_left, left_norms = _normalize_columns(left) + if numpy.any(right_norms == 0) or numpy.any(left_norms == 0): + biorthogonality_condition = numpy.inf + else: + normalized_gram = normalized_left.T @ normalized_right + biorthogonality_condition = float(numpy.linalg.cond(normalized_gram)) + powers = numpy.asarray([ + waveguide_2d.inner_product(e_field, h_field, basis.overlap_dxes, conj_h=True) + for e_field, h_field in basis.modes + ]) + full_weights = numpy.tile(numpy.outer( + 0.5 * (numpy.real(basis.overlap_dxes[0][0]) + numpy.real(basis.overlap_dxes[1][0])), + 0.5 * (numpy.real(basis.overlap_dxes[0][1]) + numpy.real(basis.overlap_dxes[1][1])), + ).ravel(), 4) + pml_weights = full_weights * numpy.tile(~basis.physical_mask.ravel(), 4) + pml_participation = numpy.sum(pml_weights[:, None] * numpy.abs(right) ** 2, axis=0) / numpy.sum( + full_weights[:, None] * numpy.abs(right) ** 2, + axis = 0, + ) + diagnostics = dataclasses.replace( + basis.diagnostics, + right_residuals = _block_relative_residuals(operator, basis.wavenumber_matrix, right), + left_residuals = _block_relative_residuals(operator.T, basis.wavenumber_matrix.T, left), + powers = powers, + pml_participation = numpy.asarray(pml_participation, dtype=float), + biorthogonality_error = float(numpy.max(numpy.abs(gram - numpy.eye(len(basis))))), + biorthogonality_condition = biorthogonality_condition, + ) + return dataclasses.replace(basis, diagnostics=diagnostics) + + +def _transform_basis( + basis: ComplexModeBasis, + transform: NDArray[numpy.complex128], + *, + spectrum_order: NDArray[numpy.int_] | None = None, + ) -> ComplexModeBasis: + """Apply a modal basis transform and its contragredient dual transform.""" + transform = numpy.asarray(transform, dtype=complex) + if transform.shape != (len(basis), len(basis)) or not numpy.isfinite(transform).all(): + raise ValueError('basis transform must be a finite square mode matrix') + inverse = numpy.linalg.inv(transform) + diagnostics = basis.diagnostics + wavenumbers = basis.wavenumbers + if spectrum_order is not None: + wavenumbers = wavenumbers[spectrum_order] + diagnostics = dataclasses.replace( + diagnostics, + target_indices = diagnostics.target_indices[spectrum_order], + target_distances = diagnostics.target_distances[spectrum_order], + ) + transformed = dataclasses.replace( + basis, + wavenumbers = wavenumbers, + wavenumber_matrix = inverse @ basis.wavenumber_matrix @ transform, + electric_fields = transform.T @ basis.electric_fields, + magnetic_fields = transform.T @ basis.magnetic_fields, + dual_electric_fields = inverse @ basis.dual_electric_fields, + dual_magnetic_fields = inverse @ basis.dual_magnetic_fields, + diagnostics = diagnostics, + ) + return _refresh_diagnostics(transformed) + + +def _scaled_basis(basis: ComplexModeBasis, scales: NDArray[numpy.complex128]) -> ComplexModeBasis: + """Scale primal basis vectors and reciprocally scale their Lorentz duals.""" + return _transform_basis(basis, numpy.diag(scales)) + + +def power_normalize_modes( + basis: ComplexModeBasis, + indices: Sequence[int], + *, + power_atol: float = 1e-12, + power_rtol: float = 1e-6, + ) -> ComplexModeBasis: + """Normalize selected positive-forward-power modes to unit power. + + Primal fields are scaled by the inverse square root of their real power and + Lorentz duals are scaled reciprocally, preserving biorthogonality and the + represented beta operator. + + Args: + basis: Modal basis to copy and rescale. + indices: Basis-vector indices to normalize. + power_atol: Absolute threshold for positive power and negligible + imaginary power. + power_rtol: Relative threshold for negligible imaginary power. + + Returns: + A transformed basis with refreshed diagnostics. + + Raises: + ValueError: If tolerances are invalid or a selected mode does not carry + positive, nearly real forward power. + """ + if not numpy.isfinite(power_atol) or power_atol <= 0: + raise ValueError('power_atol must be finite and positive') + if not numpy.isfinite(power_rtol) or power_rtol <= 0: + raise ValueError('power_rtol must be finite and positive') + scales = numpy.ones(len(basis), dtype=complex) + for index in indices: + power = basis.diagnostics.powers[index] + if power.real <= power_atol or abs(power.imag) > max(power_atol, power_rtol * abs(power.real)): + raise ValueError(f'mode {index} does not have positive nearly-real forward power') + scales[index] = 1 / numpy.sqrt(power.real) + return _scaled_basis(basis, scales) + + +def _lorentz_overlap( + dual_mode: tuple[vcfdfield2, vcfdfield2], + mode: tuple[vcfdfield2, vcfdfield2], + dxes: dx_lists2_t, + ) -> complex: + """Evaluate a dual/primal symmetric unconjugated Lorentz overlap.""" + return 0.5 * ( + waveguide_2d.inner_product(dual_mode[0], mode[1], dxes, conj_h=False) + + waveguide_2d.inner_product(mode[0], dual_mode[1], dxes, conj_h=False) + ) + + +def _reorder_basis(basis: ComplexModeBasis, order: NDArray[numpy.int_]) -> ComplexModeBasis: + """Permute a basis and its associated spectral labels.""" + transform = numpy.eye(len(basis), dtype=complex)[:, order] + return _transform_basis(basis, transform, spectrum_order=order) + + +def _degenerate_clusters( + previous: NDArray[numpy.complex128], + current: NDArray[numpy.complex128], + rtol: float, + ) -> list[list[int]]: + """Find connected mode groups degenerate in two adjacent spectra.""" + count = previous.size + remaining = set(range(count)) + clusters: list[list[int]] = [] + while remaining: + seed = remaining.pop() + cluster = {seed} + changed = True + while changed: + changed = False + for candidate in tuple(remaining): + if any( + abs(previous[candidate] - previous[item]) <= rtol * max(1.0, abs(previous[candidate]), abs(previous[item])) + and abs(current[candidate] - current[item]) <= rtol * max(1.0, abs(current[candidate]), abs(current[item])) + for item in cluster + ): + remaining.remove(candidate) + cluster.add(candidate) + changed = True + clusters.append(sorted(cluster)) + return clusters + + +def track_mode_path( + bases: Sequence[ComplexModeBasis], + *, + degeneracy_rtol: float = 1e-6, + minimum_overlap: float = 1e-3, + ) -> tuple[ComplexModeBasis, ...]: + """Assign, align, and parallel-transport modal bases along a path. + + Adjacent modes are assigned by Lorentz overlap, with beta distance used as + a deterministic tie-breaker. Nearly degenerate groups are aligned by a + unitary Procrustes rotation, so their propagation operator may become + non-diagonal in the transported basis. + + Args: + bases: Equal-sized modal bases sampled in path order. + degeneracy_rtol: Relative spectral tolerance used to group modes for + subspace rotation. + minimum_overlap: Minimum accepted adjacent diagonal overlap after + assignment and subspace alignment. + + Returns: + Basis sequence in tracked order with aligned phases, contragredient + duals, transformed beta matrices, and refreshed diagnostics. + + Raises: + ValueError: If the sequence, tolerances, field shapes, or overlap grids + are incompatible. + RuntimeError: If an assigned mode has insufficient adjacent overlap. + """ + if not bases: + raise ValueError('bases must not be empty') + if not numpy.isfinite(degeneracy_rtol) or degeneracy_rtol <= 0: + raise ValueError('degeneracy_rtol must be finite and positive') + if not numpy.isfinite(minimum_overlap) or minimum_overlap <= 0: + raise ValueError('minimum_overlap must be finite and positive') + count = len(bases[0]) + if count == 0 or any(len(basis) != count for basis in bases): + raise ValueError('all bases must contain the same nonzero mode count') + reference = bases[0] + for index, basis in enumerate(bases[1:], start=1): + if basis.electric_fields.shape != reference.electric_fields.shape or basis.magnetic_fields.shape != reference.magnetic_fields.shape: + raise ValueError(f'bases[{index}] modal field shapes do not match bases[0]') + for family in range(2): + for axis in range(2): + left_grid = numpy.asarray(reference.overlap_dxes[family][axis]) + right_grid = numpy.asarray(basis.overlap_dxes[family][axis]) + if left_grid.shape != right_grid.shape or not numpy.allclose(left_grid, right_grid, rtol=1e-12, atol=0): + raise ValueError('all bases must use identical overlap grids; resample fields before tracking') + + tracked = [bases[0]] + for current_unordered in bases[1:]: + previous = tracked[-1] + overlap = numpy.asarray([ + [_lorentz_overlap(dual_mode, mode, previous.overlap_dxes) for mode in current_unordered.modes] + for dual_mode in previous.dual_modes + ]) + beta_scale = numpy.maximum(1.0, numpy.maximum( + numpy.abs(previous.wavenumbers[:, None]), + numpy.abs(current_unordered.wavenumbers[None, :]), + )) + beta_distance = numpy.abs(previous.wavenumbers[:, None] - current_unordered.wavenumbers[None, :]) / beta_scale + cost = -numpy.abs(overlap) + 1e-9 * beta_distance + rows, columns = linear_sum_assignment(cost) + order = columns[numpy.argsort(rows)] + permutation = numpy.eye(count, dtype=complex)[:, order] + reordered_e = current_unordered.electric_fields[order] + reordered_h = current_unordered.magnetic_fields[order] + reordered_wavenumbers = current_unordered.wavenumbers[order] + + cluster_transform = numpy.eye(count, dtype=complex) + for cluster in _degenerate_clusters(previous.wavenumbers, reordered_wavenumbers, degeneracy_rtol): + if len(cluster) == 1: + continue + indices = numpy.asarray(cluster) + cross = numpy.asarray([ + [_lorentz_overlap(previous.dual_modes[ii], (reordered_e[jj], reordered_h[jj]), previous.overlap_dxes) for jj in indices] + for ii in indices + ]) + left_u, _singular_values, right_vh = numpy.linalg.svd(cross) + transform = right_vh.conj().T @ left_u.conj().T + cluster_transform[numpy.ix_(indices, indices)] = transform + + aligned_e = cluster_transform.T @ reordered_e + aligned_h = cluster_transform.T @ reordered_h + diagonal_overlap = numpy.asarray([ + _lorentz_overlap(previous.dual_modes[index], (aligned_e[index], aligned_h[index]), previous.overlap_dxes) + for index in range(count) + ]) + if numpy.any(numpy.abs(diagonal_overlap) < minimum_overlap): + bad = numpy.flatnonzero(numpy.abs(diagonal_overlap) < minimum_overlap) + raise RuntimeError(f'could not track modes {bad.tolist()}: adjacent overlap is too small') + phase_scales = numpy.exp(-1j * numpy.angle(diagonal_overlap)) + transform = permutation @ cluster_transform @ numpy.diag(phase_scales) + tracked.append(_transform_basis(current_unordered, transform, spectrum_order=order)) + + return tuple(tracked) diff --git a/meanas/test/test_eigensolvers_numerics.py b/meanas/test/test_eigensolvers_numerics.py index cf037a6..88de789 100644 --- a/meanas/test/test_eigensolvers_numerics.py +++ b/meanas/test/test_eigensolvers_numerics.py @@ -1,10 +1,11 @@ import numpy from numpy.linalg import norm import pytest +from scipy import sparse from ._solver_cases import diagonal_eigen_case from .utils import assert_close -from ..eigensolvers import power_iteration, rayleigh_quotient_iteration, signed_eigensolve +from ..eigensolvers import power_iteration, rayleigh_quotient_iteration, scipy_eigs_near, signed_eigensolve def test_rayleigh_quotient_iteration_with_linear_operator() -> None: @@ -99,3 +100,27 @@ def test_signed_eigensolve_returns_largest_positive_modes() -> None: assert_close(eigvals, [3.0, 5.0], atol=1e-6) assert eigvecs.shape == (4, 2) + + +def test_scipy_eigs_near_returns_paired_bilinear_left_vectors() -> None: + operator = sparse.csc_array(numpy.diag([1.0 + 0.1j, 2.0 - 0.2j, 4.0 + 0.3j, 7.0])) + + solution = scipy_eigs_near(operator, target=3.1, how_many=2) + + assert_close(numpy.sort_complex(solution.eigenvalues), numpy.sort_complex([2.0 - 0.2j, 4.0 + 0.3j])) + gram = solution.left_eigenvectors.T @ solution.right_eigenvectors + assert numpy.linalg.cond(gram) < 2 + for index, eigenvalue in enumerate(solution.eigenvalues): + right = solution.right_eigenvectors[:, index] + left = solution.left_eigenvectors[:, index] + assert norm(operator @ right - eigenvalue * right) < 1e-10 + assert norm(operator.T @ left - eigenvalue * left) < 1e-10 + + +@pytest.mark.parametrize('parameter', ['tolerance', 'pairing_rtol']) +@pytest.mark.parametrize('invalid', [numpy.nan, numpy.inf, -numpy.inf]) +def test_scipy_eigs_near_rejects_nonfinite_tolerances(parameter: str, invalid: float) -> None: + operator = sparse.csc_array(numpy.diag([1.0, 2.0, 3.0, 4.0])) + + with pytest.raises(ValueError, match=rf'{parameter} must be finite and positive'): + scipy_eigs_near(operator, target=2.1, how_many=1, **{parameter: invalid}) diff --git a/meanas/test/test_waveguide_2d_complex.py b/meanas/test/test_waveguide_2d_complex.py new file mode 100644 index 0000000..fc1fb1c --- /dev/null +++ b/meanas/test/test_waveguide_2d_complex.py @@ -0,0 +1,694 @@ +import dataclasses +import numpy +import pytest +from numpy.testing import assert_allclose +from scipy import sparse +from scipy import linalg + +from ..eigensolvers import TargetedEigensolution, scipy_eigs_near +from ..fdmath import vec +from ..fdfd import eme, scpml, waveguide_2d, waveguide_2d_complex, waveguide_3d + + +def _grid(shape: tuple[int, int], spacing: float = 0.3) -> list[list[numpy.ndarray]]: + return [[numpy.full(shape[0], spacing), numpy.full(shape[1], spacing)] for _ in range(2)] + + +def _waveguide_case( + shape: tuple[int, int] = (11, 11), + *, + pml_cells: int = 0, + ) -> tuple[float, list[list[numpy.ndarray]], numpy.ndarray]: + omega = 2 * numpy.pi / 1.55 + epsilon_cladding = 1.45 ** 2 + epsilon = numpy.full((3, *shape), epsilon_cladding) + cx, cy = (size // 2 for size in shape) + epsilon[:, cx - 1:cx + 2, cy - 1:cy + 2] = 2.0 ** 2 + dxes = _grid(shape) + if pml_cells: + for axis in (0, 1): + for polarity in (-1, +1): + dxes = scpml.stretch_with_scpml( + dxes, + axis=axis, + polarity=polarity, + omega=omega, + epsilon_effective=epsilon_cladding, + thickness=pml_cells, + ) + return omega, dxes, vec(epsilon) + + +def test_operator_beta_squares_to_existing_e_and_h_operators() -> None: + rng = numpy.random.default_rng(4) + shape = (4, 5) + omega = 1.7 + 0.03j + dxes = [[ + 0.7 + rng.random(size) + 0.02j * rng.random(size) + for size in shape + ] for _ in range(2)] + epsilon = vec(1.5 + rng.random((3, *shape))) + mu = vec(0.9 + 0.2 * rng.random((3, *shape))) + + operator = waveguide_2d_complex.operator_beta(omega, dxes, epsilon, mu) + n_state = operator.shape[0] // 2 + squared = operator @ operator + electric = squared[:n_state, :n_state] + magnetic = squared[n_state:, n_state:] + + assert_allclose(electric.toarray(), waveguide_2d.operator_e(omega, dxes, epsilon, mu).toarray(), atol=3e-12) + assert_allclose(magnetic.toarray(), waveguide_2d.operator_h(omega, dxes, epsilon, mu).toarray(), atol=3e-12) + assert sparse.linalg.norm(squared[:n_state, n_state:]) == 0 + assert sparse.linalg.norm(squared[n_state:, :n_state]) == 0 + + +def test_complex_solver_matches_guided_beta_and_builds_lorentz_duals() -> None: + omega, dxes, epsilon = _waveguide_case() + _e_xy, expected_beta = waveguide_2d.solve_mode(0, omega=omega, dxes=dxes, epsilon=epsilon) + + basis = waveguide_2d_complex.solve_modes_near( + omega, + dxes, + epsilon, + target_wavenumbers=[0.98 * expected_beta], + modes_per_target=1, + candidate_margin=2, + ) + + assert_allclose(basis.wavenumbers[0], expected_beta, rtol=1e-4, atol=1e-10) + assert basis.diagnostics.right_residuals[0] < 1e-10 + assert basis.diagnostics.left_residuals[0] < 1e-10 + overlap = 0.5 * ( + waveguide_2d.inner_product(basis.dual_electric_fields[0], basis.magnetic_fields[0], dxes) + + waveguide_2d.inner_product(basis.electric_fields[0], basis.dual_magnetic_fields[0], dxes) + ) + assert_allclose(overlap, 1.0, atol=1e-10) + + +def test_complex_solver_returns_passive_pml_modes_and_diagnostics() -> None: + omega, dxes, epsilon = _waveguide_case(shape=(15, 15), pml_cells=3) + + basis = waveguide_2d_complex.solve_modes_near( + omega, + dxes, + epsilon, + target_wavenumbers=[1.45 * omega], + modes_per_target=3, + candidate_margin=2, + ) + + assert numpy.all(numpy.imag(basis.wavenumbers) < 1e-8) + assert numpy.all((basis.diagnostics.pml_participation >= 0) & (basis.diagnostics.pml_participation <= 1)) + assert basis.diagnostics.biorthogonality_error < 1e-10 + + +def test_mode_tracking_restores_permutation_phase_and_eme_adapter() -> None: + omega, dxes, epsilon = _waveguide_case(shape=(15, 15), pml_cells=3) + basis = waveguide_2d_complex.solve_modes_near( + omega, + dxes, + epsilon, + target_wavenumbers=[1.45 * omega], + modes_per_target=2, + candidate_margin=2, + ) + permuted = waveguide_2d_complex._reorder_basis(basis, numpy.array([1, 0])) + permuted = waveguide_2d_complex._scaled_basis(permuted, numpy.exp(1j * numpy.array([0.7, -0.4]))) + + tracked = waveguide_2d_complex.track_mode_path([basis, permuted]) + assert_allclose(tracked[1].wavenumbers, basis.wavenumbers) + for index in range(len(basis)): + overlap = waveguide_2d_complex._lorentz_overlap( + tracked[0].dual_modes[index], + tracked[1].modes[index], + tracked[0].overlap_dxes, + ) + assert overlap.real > 0 + assert abs(overlap.imag) < 1e-10 + + sections = eme.mode_sections_from_bases([0.0, 2.0], tracked) + assert len(sections) == 2 + assert len(sections[0].modes) == len(basis) + assert sections[0].dual_modes is not None + assert sections[0].backward_modes is None + + +def test_power_normalization_preserves_dual_overlap() -> None: + omega, dxes, epsilon = _waveguide_case() + _e_xy, beta = waveguide_2d.solve_mode(0, omega=omega, dxes=dxes, epsilon=epsilon) + basis = waveguide_2d_complex.solve_modes_near( + omega, + dxes, + epsilon, + target_wavenumbers=[0.98 * beta], + modes_per_target=1, + candidate_margin=2, + ) + + normalized = waveguide_2d_complex.power_normalize_modes(basis, [0]) + + assert_allclose(normalized.diagnostics.powers[0], 1.0, atol=1e-10) + overlap = waveguide_2d_complex._lorentz_overlap(normalized.dual_modes[0], normalized.modes[0], dxes) + assert_allclose(overlap, 1.0, atol=1e-10) + for parameter in ('power_atol', 'power_rtol'): + for invalid in (numpy.nan, numpy.inf, -numpy.inf): + with pytest.raises(ValueError, match=rf'{parameter} must be finite and positive'): + waveguide_2d_complex.power_normalize_modes(basis, [0], **{parameter: invalid}) + + +def test_overlapping_targets_merge_instead_of_backfilling_farther_modes() -> None: + values = numpy.arange(1, 5, dtype=complex) + vectors = numpy.eye(4, dtype=complex) + solution = TargetedEigensolution(values, vectors, vectors) + + selected, _right, _left, target_indices = waveguide_2d_complex._select_unique_modes( + [solution, solution], + numpy.asarray([1.1, 1.2]), + 2, + eigenvalue_rtol=1e-8, + subspace_tol=1e-8, + state_weights=numpy.ones(4), + ) + + assert_allclose(selected, [1, 2]) + assert_allclose(target_indices, [0, 1]) + + +def test_rotated_degenerate_candidate_subspaces_are_deduplicated_by_rank() -> None: + values = numpy.asarray([2.0, 2.0]) + vectors = numpy.eye(4, 2, dtype=complex) + rotation = numpy.asarray([[1, 1j], [1j, 1]], dtype=complex) / numpy.sqrt(2) + first = TargetedEigensolution(values, vectors, vectors) + second = TargetedEigensolution(values, vectors @ rotation, vectors @ rotation) + + selected, right, left, _targets = waveguide_2d_complex._select_unique_modes( + [first, second], + numpy.asarray([1.9, 2.1]), + 2, + eigenvalue_rtol=1e-8, + subspace_tol=1e-8, + state_weights=numpy.ones(4), + ) + + assert selected.size == 2 + assert numpy.linalg.matrix_rank(right) == 2 + assert numpy.linalg.matrix_rank(left) == 2 + + +def test_independently_cross_paired_degenerate_subspaces_are_accepted() -> None: + e1 = numpy.asarray([1.0, 0.0], dtype=complex) + e2 = numpy.asarray([0.0, 1.0], dtype=complex) + first = TargetedEigensolution(numpy.asarray([2.0]), numpy.column_stack([e1]), numpy.column_stack([e1])) + second = TargetedEigensolution( + numpy.asarray([2.0, 2.0]), + numpy.column_stack([e2, e1]), + numpy.column_stack([e1, e2]), + ) + + selected, right, left, _targets = waveguide_2d_complex._select_unique_modes( + [first, second], numpy.asarray([1.9, 2.1]), 1, + eigenvalue_rtol=1e-8, subspace_tol=1e-8, state_weights=numpy.ones(2), + ) + + assert selected.size == 2 + assert numpy.linalg.matrix_rank(right) == 2 + assert numpy.linalg.matrix_rank(left) == 2 + + +def test_single_target_cross_paired_degenerate_candidates_form_dual_basis() -> None: + e1 = numpy.asarray([1.0, 0.0], dtype=complex) + e2 = numpy.asarray([0.0, 1.0], dtype=complex) + solution = TargetedEigensolution( + numpy.asarray([2.0, 2.0]), + numpy.column_stack([e1, e2]), + numpy.column_stack([e2, e1]), + ) + + selected, right, left, _targets = waveguide_2d_complex._select_unique_modes( + [solution], numpy.asarray([2.1]), 1, + eigenvalue_rtol=1e-8, subspace_tol=1e-8, state_weights=numpy.ones(2), + ) + + assert selected.size == 1 + assert right.shape == (2, 1) + assert left.shape == (2, 1) + assert abs((left.T @ right)[0, 0]) > 0.5 + + +def test_reserve_candidate_repairs_higher_rank_cross_gram() -> None: + vectors = numpy.eye(3, dtype=complex) + solution = TargetedEigensolution( + numpy.full(3, 2.0), + vectors, + vectors[:, [1, 2, 0]], + ) + + selected, right, left, _targets = waveguide_2d_complex._select_unique_modes( + [solution], numpy.asarray([2.1]), 2, + eigenvalue_rtol=1e-8, subspace_tol=1e-8, state_weights=numpy.ones(3), + ) + + assert selected.size == 2 + assert numpy.linalg.matrix_rank(left.T @ right) == 2 + + +def test_cluster_with_insufficient_cross_rank_is_rejected() -> None: + e1 = numpy.asarray([1.0, 0.0], dtype=complex) + e2 = numpy.asarray([0.0, 1.0], dtype=complex) + first = TargetedEigensolution(numpy.asarray([2.0]), numpy.column_stack([e1]), numpy.column_stack([e1])) + second = TargetedEigensolution(numpy.asarray([2.0]), numpy.column_stack([e2]), numpy.column_stack([e1])) + + with pytest.raises(RuntimeError, match=r'required cross rank 2, available 1'): + waveguide_2d_complex._select_unique_modes( + [first, second], numpy.asarray([1.9, 2.1]), 1, + eigenvalue_rtol=1e-8, subspace_tol=1e-8, state_weights=numpy.ones(2), + ) + + +def test_single_target_does_not_expand_to_full_degenerate_candidate_margin() -> None: + vectors = numpy.eye(5, dtype=complex) + solution = TargetedEigensolution(numpy.full(5, 2.0), vectors, vectors) + + selected, right, left, _targets = waveguide_2d_complex._select_unique_modes( + [solution], numpy.asarray([2.1]), 1, + eigenvalue_rtol=1e-8, subspace_tol=1e-8, state_weights=numpy.ones(5), + ) + + assert selected.size == 1 + assert right.shape == (5, 1) + assert left.shape == (5, 1) + + +def test_mode_selection_is_invariant_to_candidate_scaling() -> None: + values = numpy.asarray([1.0, 2.0], dtype=complex) + vectors = numpy.eye(2, dtype=complex) + reference = TargetedEigensolution(values, vectors, vectors) + scaled = TargetedEigensolution( + values, + vectors * numpy.asarray([1e-12, 1e8]), + vectors * numpy.asarray([1e9j, -1e-11]), + ) + + outputs = [ + waveguide_2d_complex._select_unique_modes( + [solution], numpy.asarray([1.5]), 2, + eigenvalue_rtol=1e-8, subspace_tol=1e-7, state_weights=numpy.ones(2), + ) + for solution in (reference, scaled) + ] + + assert_allclose(outputs[0][0], outputs[1][0]) + for _values, right, left, _targets in outputs: + assert numpy.linalg.matrix_rank(left.T @ right) == 2 + + +def test_mode_selection_is_invariant_to_uniform_weight_scaling() -> None: + values = numpy.asarray([1.0, 2.0], dtype=complex) + vectors = numpy.eye(2, dtype=complex) + solution = TargetedEigensolution(values, vectors, vectors) + + outputs = [ + waveguide_2d_complex._select_unique_modes( + [solution], numpy.asarray([1.5]), 2, + eigenvalue_rtol=1e-8, subspace_tol=1e-7, + state_weights=numpy.full(2, weight_scale), + ) + for weight_scale in (1.0, 1e8) + ] + + assert_allclose(outputs[0][0], outputs[1][0]) + assert_allclose(outputs[0][1], outputs[1][1]) + assert_allclose(outputs[0][2], outputs[1][2]) + + +def test_complex_solver_is_invariant_to_backend_eigenvector_scaling() -> None: + omega, dxes, epsilon = _waveguide_case() + eigensolver_opts = {'v0': numpy.linspace(1.0, 2.0, 4 * 11 * 11, dtype=complex)} + cached_solution = None + + def cached_backend(operator, target, count, **options): + nonlocal cached_solution + if cached_solution is None: + cached_solution = scipy_eigs_near(operator, target, count, **options) + return cached_solution + + reference = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=2, candidate_margin=2, + eigensolver=cached_backend, eigensolver_opts=eigensolver_opts, + ) + + def rescaled_backend(operator, target, count, **options): + solution = cached_backend(operator, target, count, **options) + order = numpy.argsort(numpy.abs(solution.eigenvalues - target)) + right_scales = numpy.ones(count) + right_scales[order[:2]] = [1e-16, 1e16] + return TargetedEigensolution( + solution.eigenvalues, + solution.right_eigenvectors * right_scales, + solution.left_eigenvectors, + ) + + scaled = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=2, candidate_margin=2, + eigensolver=rescaled_backend, eigensolver_opts=eigensolver_opts, + ) + + assert_allclose( + numpy.sort_complex(scaled.wavenumbers), + numpy.sort_complex(reference.wavenumbers), + atol=1e-10, + rtol=1e-10, + ) + assert numpy.max(scaled.diagnostics.right_residuals) < 1e-8 + assert numpy.max(scaled.diagnostics.left_residuals) < 1e-8 + assert scaled.diagnostics.biorthogonality_error < 1e-9 + assert_allclose( + scaled.diagnostics.biorthogonality_condition, + reference.diagnostics.biorthogonality_condition, + rtol=1e-10, + atol=1e-12, + ) + + +def test_refreshed_biorthogonality_condition_is_invariant_to_reciprocal_scaling() -> None: + omega, dxes, epsilon = _waveguide_case() + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=2, candidate_margin=2, + ) + + scaled = waveguide_2d_complex._scaled_basis(basis, numpy.asarray([1e-4, 1e4])) + + assert_allclose( + scaled.diagnostics.biorthogonality_condition, + basis.diagnostics.biorthogonality_condition, + rtol=1e-10, + atol=1e-12, + ) + assert scaled.diagnostics.biorthogonality_error < 1e-7 + + +@pytest.mark.parametrize('parameter', [ + 'residual_tol', + 'duplicate_eigenvalue_rtol', + 'duplicate_subspace_tol', +]) +@pytest.mark.parametrize('invalid', [numpy.nan, numpy.inf, -numpy.inf]) +def test_complex_solver_rejects_nonfinite_tolerances(parameter: str, invalid: float) -> None: + omega, dxes, epsilon = _waveguide_case() + + with pytest.raises(ValueError, match=rf'{parameter} must be finite and positive'): + waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=1, + **{parameter: invalid}, + ) + + +def test_basis_rotation_carries_full_wavenumber_matrix_and_residuals() -> None: + omega, dxes, epsilon = _waveguide_case(shape=(13, 13), pml_cells=2) + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=2, candidate_margin=2, + ) + rotation = numpy.asarray([[1, 1j], [1j, 1]], dtype=complex) / numpy.sqrt(2) + + rotated = waveguide_2d_complex._transform_basis(basis, rotation) + expected = numpy.linalg.inv(rotation) @ numpy.diag(basis.wavenumbers) @ rotation + + assert_allclose(rotated.wavenumber_matrix, expected, atol=1e-11) + assert_allclose(numpy.sort_complex(numpy.linalg.eigvals(rotated.wavenumber_matrix)), numpy.sort_complex(basis.wavenumbers), atol=1e-10) + assert numpy.max(rotated.diagnostics.right_residuals) < 1e-8 + assert numpy.max(rotated.diagnostics.left_residuals) < 1e-8 + assert rotated.diagnostics.biorthogonality_error < 1e-9 + + +def test_eme_propagates_constant_mixed_basis_with_full_matrix() -> None: + omega, dxes, epsilon = _waveguide_case(shape=(13, 13), pml_cells=2) + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=2, candidate_margin=2, + ) + rotation = numpy.asarray([[1, 1j], [1j, 1]], dtype=complex) / numpy.sqrt(2) + rotated = waveguide_2d_complex._transform_basis(basis, rotation) + length = 0.7 + sections = eme.mode_sections_from_bases([0.0, length], [rotated, rotated]) + + transfer = eme.get_taper_abcd(sections, dxes=rotated.overlap_dxes).toarray() + expected = linalg.block_diag( + linalg.expm(-1j * rotated.wavenumber_matrix * length), + linalg.expm(+1j * rotated.wavenumber_matrix * length), + ) + assert_allclose(transfer, expected, atol=1e-9, rtol=1e-9) + + +def test_eme_adapter_rejects_subset_that_cuts_coupled_block() -> None: + omega, dxes, epsilon = _waveguide_case(shape=(13, 13), pml_cells=2) + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=2, candidate_margin=2, + ) + rotation = numpy.asarray([[1, 1j], [1j, 1]], dtype=complex) / numpy.sqrt(2) + rotated = waveguide_2d_complex._transform_basis(basis, rotation) + + with pytest.raises(ValueError, match='complete coupled block'): + eme.mode_sections_from_bases([0.0, 1.0], [rotated, rotated], mode_indices=[0]) + with pytest.raises(ValueError, match='complete coupled block'): + eme.mode_sections_from_bases( + [0.0, 1.0], [basis, basis], mode_indices=[0], backward_bases=[rotated, rotated], + ) + sections = eme.mode_sections_from_bases([0.0, 1.0], [basis, basis], mode_indices=[0]) + assert len(sections[0].modes) == 1 + for invalid_tolerance in (numpy.nan, numpy.inf, -numpy.inf): + with pytest.raises(ValueError, match='finite and positive'): + eme.mode_sections_from_bases([0.0, 1.0], [basis, basis], coupling_rtol=invalid_tolerance) + for invalid_indices in ([0.9], [True]): + with pytest.raises(ValueError, match='entries must be integers'): + eme.mode_sections_from_bases([0.0, 1.0], [basis, basis], mode_indices=invalid_indices) + + +def test_mode_tracking_rejects_different_overlap_grids() -> None: + omega, dxes, epsilon = _waveguide_case() + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=1, candidate_margin=2, + ) + changed_dxes = [[values.copy() for values in family] for family in basis.overlap_dxes] + changed_dxes[0][0][0] *= 1.01 + changed = dataclasses.replace(basis, overlap_dxes=changed_dxes) + + with pytest.raises(ValueError, match='identical overlap grids'): + waveguide_2d_complex.track_mode_path([basis, changed]) + for parameter in ('degeneracy_rtol', 'minimum_overlap'): + for invalid in (numpy.nan, numpy.inf, -numpy.inf): + with pytest.raises(ValueError, match=rf'{parameter} must be finite and positive'): + waveguide_2d_complex.track_mode_path([basis], **{parameter: invalid}) + + +def test_mode_tracking_refreshes_diagnostics_once_per_added_section(monkeypatch: pytest.MonkeyPatch) -> None: + omega, dxes, epsilon = _waveguide_case(shape=(13, 13), pml_cells=2) + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=2, candidate_margin=2, + ) + permuted = dataclasses.replace( + basis, + wavenumbers=basis.wavenumbers[::-1], + wavenumber_matrix=basis.wavenumber_matrix[::-1, ::-1], + electric_fields=basis.electric_fields[::-1], + magnetic_fields=basis.magnetic_fields[::-1], + dual_electric_fields=basis.dual_electric_fields[::-1], + dual_magnetic_fields=basis.dual_magnetic_fields[::-1], + ) + original_refresh = waveguide_2d_complex._refresh_diagnostics + refresh_count = 0 + + def counted_refresh(current_basis): + nonlocal refresh_count + refresh_count += 1 + return original_refresh(current_basis) + + monkeypatch.setattr(waveguide_2d_complex, '_refresh_diagnostics', counted_refresh) + tracked = waveguide_2d_complex.track_mode_path([basis, permuted]) + + assert refresh_count == 1 + assert_allclose(numpy.sort_complex(numpy.linalg.eigvals(tracked[1].wavenumber_matrix)), numpy.sort_complex(basis.wavenumbers), atol=1e-9) + assert tracked[1].diagnostics.biorthogonality_error < 1e-9 + + +def test_mode_tracking_procrustes_path_carries_nondiagonal_wavenumber_matrix() -> None: + omega, dxes, epsilon = _waveguide_case(shape=(13, 13), pml_cells=2) + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=2, candidate_margin=2, + ) + rotation = numpy.asarray([[1, 1j], [1j, 1]], dtype=complex) / numpy.sqrt(2) + previous = waveguide_2d_complex._transform_basis(basis, rotation) + + tracked = waveguide_2d_complex.track_mode_path([previous, basis], degeneracy_rtol=1.0) + off_diagonal = tracked[1].wavenumber_matrix - numpy.diag(numpy.diag(tracked[1].wavenumber_matrix)) + + assert numpy.linalg.norm(off_diagonal) > 1e-8 + assert_allclose(tracked[1].wavenumber_matrix, previous.wavenumber_matrix, atol=1e-9, rtol=1e-9) + assert numpy.max(tracked[1].diagnostics.right_residuals) < 1e-8 + assert numpy.max(tracked[1].diagnostics.left_residuals) < 1e-8 + assert tracked[1].diagnostics.biorthogonality_error < 1e-9 + + +def test_signed_beta_targets_return_forward_backward_pair() -> None: + omega, dxes, epsilon = _waveguide_case() + _field, beta = waveguide_2d.solve_mode(0, omega=omega, dxes=dxes, epsilon=epsilon) + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[0.99 * beta, -0.99 * beta], modes_per_target=1, candidate_margin=2, + ) + + assert basis.wavenumbers.size == 2 + assert_allclose(basis.wavenumbers[0], -basis.wavenumbers[1], rtol=1e-8, atol=1e-9) + n_points = 11 * 11 + phase = numpy.vdot(basis.electric_fields[0, :2 * n_points], basis.electric_fields[1, :2 * n_points]) + phase /= abs(phase) + assert_allclose(basis.electric_fields[0, :2 * n_points], basis.electric_fields[1, :2 * n_points] / phase, rtol=1e-7, atol=1e-8) + assert_allclose(basis.magnetic_fields[0, :2 * n_points], -basis.magnetic_fields[1, :2 * n_points] / phase, rtol=1e-7, atol=1e-8) + + +def test_guided_mode_is_stable_under_window_and_pml_refinement() -> None: + results = [] + for shape, pml_cells in ((15, 3), (19, 4)): + omega, dxes, epsilon = _waveguide_case(shape=(shape, shape), pml_cells=pml_cells) + basis = waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.8 * omega], modes_per_target=1, candidate_margin=2, + ) + results.append(basis) + + assert_allclose(results[0].wavenumbers, results[1].wavenumbers, rtol=1e-4, atol=1e-8) + assert results[1].diagnostics.pml_participation[0] < results[0].diagnostics.pml_participation[0] + + +def test_invalid_targeted_backend_result_is_rejected() -> None: + omega, dxes, epsilon = _waveguide_case() + + def invalid_backend(operator, _target, _count): + return TargetedEigensolution( + numpy.asarray([1.0]), + numpy.ones((operator.shape[0] - 1, 1)), + numpy.ones((operator.shape[0] - 1, 1)), + ) + + with pytest.raises(ValueError, match='eigenvector shapes'): + waveguide_2d_complex.solve_modes_near( + omega, dxes, epsilon, + target_wavenumbers=[1.45 * omega], modes_per_target=1, eigensolver=invalid_backend, + ) + + +@pytest.mark.complete +def test_richer_complex_basis_reduces_driven_fdfd_monitor_residual() -> None: + omega = 2 * numpy.pi / 1.55 + shape = (7, 7) + dxes_2d = _grid(shape) + left = numpy.full((3, *shape), 1.45 ** 2) + right = left.copy() + left[:, 2:5, 3:4] = 2.0 ** 2 + right[:, 2:5, 2:5] = 2.0 ** 2 + _ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports( + left, right, + omega=omega, + dxes_2d=dxes_2d, + mode_numbers=[0], + axis=2, + pml_cells=4, + pml_buffer_cells=2, + source_monitor_cells=4, + monitor_interface_cells=4, + transverse_padding_cells=2, + transverse_pml_cells=2, + pml_epsilon_effective=1.45 ** 2, + matrix_solver_opts={'atol': 1e-10, 'rtol': 1e-8}, + return_diagnostics=True, + ) + monitor_index = diagnostics['left_monitor_slice'][2].start + e_field = diagnostics['electric_fields'][0][:, :, :, monitor_index].reshape(3, -1) + h_field = diagnostics['magnetic_fields'][0][:, :, :, monitor_index].reshape(3, -1) + solve_dxes = [[diagnostics['solve_dxes'][family][axis] for axis in (0, 1)] for family in range(2)] + overlap_dxes = [[diagnostics['dxes'][family][axis] for axis in (0, 1)] for family in range(2)] + epsilon = vec(diagnostics['epsilon'][:, :, :, monitor_index]) + + def reconstruction_error(modes_per_target: int) -> float: + basis = waveguide_2d_complex.solve_modes_near( + omega, + solve_dxes, + epsilon, + overlap_dxes=overlap_dxes, + target_wavenumbers=[-1.8 * omega, -1.45 * omega, 1.45 * omega, 1.8 * omega], + modes_per_target=modes_per_target, + candidate_margin=2, + eigensolver_opts={ + 'v0': numpy.linspace(1.0, 2.0, 4 * e_field.shape[1], dtype=complex), + 'ncv': 24, + }, + ) + coefficients = numpy.asarray([ + waveguide_2d_complex._lorentz_overlap(dual_mode, (e_field.ravel(), h_field.ravel()), overlap_dxes) + for dual_mode in basis.dual_modes + ]) + reconstructed_e = coefficients @ basis.electric_fields + reconstructed_h = coefficients @ basis.magnetic_fields + mask = numpy.tile(basis.physical_mask.ravel(), 3) + residual = numpy.concatenate((e_field.ravel()[mask] - reconstructed_e[mask], h_field.ravel()[mask] - reconstructed_h[mask])) + reference = numpy.concatenate((e_field.ravel()[mask], h_field.ravel()[mask])) + return float(numpy.linalg.norm(residual) / numpy.linalg.norm(reference)) + + sparse_error = reconstruction_error(1) + rich_error = reconstruction_error(2) + assert rich_error < sparse_error + + +@pytest.mark.complete +def test_thin_core_taper_guided_loss_converges_with_basis_and_pml() -> None: + omega = 2 * numpy.pi / 1.55 + epsilon_cladding = 1.45 ** 2 + + def guided_loss(shape: int, pml_cells: int, modes_per_target: int) -> float: + overlap_dxes = _grid((shape, shape)) + solve_dxes = [[values.copy() for values in family] for family in overlap_dxes] + for axis in (0, 1): + for polarity in (-1, +1): + solve_dxes = scpml.stretch_with_scpml( + solve_dxes, axis=axis, polarity=polarity, omega=omega, + epsilon_effective=epsilon_cladding, thickness=pml_cells, + ) + center = shape // 2 + bases = [] + for edge_fill in (0.30, 0.35, 0.40): + epsilon = numpy.full((3, shape, shape), epsilon_cladding) + epsilon[:, center - 1:center + 2, center] = 2.0 ** 2 + mixed = epsilon_cladding + edge_fill * (2.0 ** 2 - epsilon_cladding) + epsilon[:, center - 1:center + 2, center - 1] = mixed + epsilon[:, center - 1:center + 2, center + 1] = mixed + basis = waveguide_2d_complex.solve_modes_near( + omega, solve_dxes, vec(epsilon), overlap_dxes=overlap_dxes, + target_wavenumbers=[1.8 * omega, 1.45 * omega], + modes_per_target=modes_per_target, candidate_margin=2, + eigensolver_opts={ + 'v0': numpy.linspace(1.0, 2.0, 4 * shape * shape, dtype=complex), + }, + ) + bases.append(waveguide_2d_complex.power_normalize_modes(basis, [0])) + tracked = waveguide_2d_complex.track_mode_path(bases) + scattering = eme.get_taper_s( + eme.mode_sections_from_bases([0.0, 5.0, 10.0], tracked), + dxes=overlap_dxes, + ) + mode_count = len(tracked[0]) + return float(1 - abs(scattering[0, 0]) ** 2 - abs(scattering[mode_count, 0]) ** 2) + + coarse_loss = guided_loss(15, 3, 1) + refined_loss = guided_loss(19, 4, 2) + assert numpy.isfinite([coarse_loss, refined_loss]).all() + assert -1e-6 <= coarse_loss <= 1 + assert -1e-6 <= refined_loss <= 1 + assert abs(coarse_loss - refined_loss) <= max(2e-5, 0.05 * abs(refined_loss))