[fdfd.waveguide_2d_complex] add 2d complex modesolver
This commit is contained in:
Forgejo Actions
parent
e851e25568
commit
0e89fe0460
9 changed files with 2021 additions and 14 deletions
|
|
@ -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.
|
||||
|
||||
|
|
|
|||
|
|
@ -2,6 +2,8 @@
|
|||
|
||||
::: meanas.fdfd.waveguide_2d
|
||||
|
||||
::: meanas.fdfd.waveguide_2d_complex
|
||||
|
||||
::: meanas.fdfd.waveguide_3d
|
||||
|
||||
::: meanas.fdfd.waveguide_cyl
|
||||
|
|
|
|||
90
examples/eme_taper_pml_modes.py
Normal file
90
examples/eme_taper_pml_modes.py
Normal file
|
|
@ -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()
|
||||
|
|
@ -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]
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
881
meanas/fdfd/waveguide_2d_complex.py
Normal file
881
meanas/fdfd/waveguide_2d_complex.py
Normal file
|
|
@ -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)
|
||||
|
|
@ -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})
|
||||
|
|
|
|||
694
meanas/test/test_waveguide_2d_complex.py
Normal file
694
meanas/test/test_waveguide_2d_complex.py
Normal file
|
|
@ -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))
|
||||
Loading…
Reference in a new issue