jan/meanas
Watch
1
0
Fork
You've already forked meanas
 0
Code Issues Pull requests Releases Wiki Activity

[fdfd.waveguide_3d] add FDFD-based S-matrix solve

Browse source
This commit is contained in:
Forgejo Actions 2026-07-14 00:21:12 -07:00
commit 37ed9d10bc
2 changed files with 762 additions and 5 deletions
Download patch file Download diff file Expand all files Collapse all files

560
meanas/fdfd/waveguide_3d.py
View file

@ -10,6 +10,8 @@ The intended workflow is:
2. Solve the corresponding 2D mode problem with `solve_mode(...)`.
3. Turn that mode into a one-sided source with `compute_source(...)`.
4. Build an overlap window with `compute_overlap_e(...)` for port readout.
5. For two-port discontinuities, `get_s_fdfd(...)` automates the source,
solve, and monitor-plane modal decomposition loop.
`polarity` is part of the public convention throughout this module:
@ -21,13 +23,14 @@ overlap window and how the expanded field is phased.
"""
from typing import Any, TypedDict, cast
import warnings
from collections.abc import Sequence
import operator
from collections.abc import Callable, Sequence
import numpy
from numpy.typing import NDArray
from numpy import complexfloating
from numpy import complexfloating, floating
from ..fdmath import vec, unvec, dx_lists_t, cfdfield_t, fdfield, cfdfield
from . import operators, waveguide_2d
from ..fdmath import vec, unvec, dx_lists2_t, dx_lists_t, cfdfield_t, fdfield, cfdfield
from . import functional, operators, scpml, solvers, waveguide_2d
class Waveguide3DMode(TypedDict):
@ -341,3 +344,552 @@ def expand_e(
E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in]
return cfdfield_t(E_expanded)
def _slice_start(name: str, slices: Sequence[slice], axis: int, shape: Sequence[int]) -> int:
if len(slices) != 3:
raise ValueError(f'{name} must contain exactly three slices')
axis_slice = slices[axis]
if axis_slice.start is None or axis_slice.stop is None:
raise ValueError(f'{name}[axis] must have explicit start and stop')
if axis_slice.step not in (None, 1):
raise ValueError(f'{name}[axis] must not use a step')
if axis_slice.stop - axis_slice.start != 1:
raise ValueError(f'{name}[axis] must select exactly one cell')
if axis_slice.start < 0 or axis_slice.stop > shape[axis]:
raise ValueError(f'{name}[axis] lies outside the simulation domain')
return axis_slice.start
def _as_slice_tuple(name: str, slices: Sequence[slice], axis: int, shape: Sequence[int]) -> tuple[slice, slice, slice]:
_slice_start(name, slices, axis, shape)
return (slices[0], slices[1], slices[2])
def _nonbool_index(value: object) -> int:
if isinstance(value, bool | numpy.bool_):
raise TypeError
return operator.index(value)
def _as_mode_numbers(name: str, mode_numbers: Sequence[int]) -> tuple[int, ...]:
modes: list[int] = []
for mode in mode_numbers:
try:
mode_i = _nonbool_index(mode)
except TypeError as err:
raise ValueError(f'{name} entries must be integers') from err
if mode_i < 0:
raise ValueError('mode numbers must be non-negative')
modes.append(mode_i)
if not modes:
raise ValueError('each FDFD port must contain at least one mode number')
return tuple(modes)
def _validate_fdfd_port_specs(
left_mode_numbers: Sequence[int],
right_mode_numbers: Sequence[int],
left_source_slice: Sequence[slice],
left_monitor_slice: Sequence[slice],
right_source_slice: Sequence[slice],
right_monitor_slice: Sequence[slice],
axis: int,
shape: Sequence[int],
) -> tuple[
tuple[int, ...],
tuple[int, ...],
tuple[slice, slice, slice],
tuple[slice, slice, slice],
tuple[slice, slice, slice],
tuple[slice, slice, slice],
]:
if axis not in (0, 1, 2):
raise ValueError('axis must be 0, 1, or 2')
left_modes_i = _as_mode_numbers('left_mode_numbers', left_mode_numbers)
right_modes_i = _as_mode_numbers('right_mode_numbers', right_mode_numbers)
left_source = _as_slice_tuple('left_source_slice', left_source_slice, axis, shape)
left_monitor = _as_slice_tuple('left_monitor_slice', left_monitor_slice, axis, shape)
right_source = _as_slice_tuple('right_source_slice', right_source_slice, axis, shape)
right_monitor = _as_slice_tuple('right_monitor_slice', right_monitor_slice, axis, shape)
left_source_i = _slice_start('left_source_slice', left_source, axis, shape)
left_monitor_i = _slice_start('left_monitor_slice', left_monitor, axis, shape)
right_monitor_i = _slice_start('right_monitor_slice', right_monitor, axis, shape)
right_source_i = _slice_start('right_source_slice', right_source, axis, shape)
if not (left_source_i < left_monitor_i < right_monitor_i < right_source_i):
raise ValueError('port slices must be ordered as left source, left monitor, right monitor, right source')
return left_modes_i, right_modes_i, left_source, left_monitor, right_source, right_monitor
def _port_decomposition(
incoming_modes: Sequence[Waveguide3DMode],
outgoing_modes: Sequence[Waveguide3DMode],
*,
dxes: dx_lists_t,
axis: int,
incoming_polarity: int,
outgoing_polarity: int,
slices: Sequence[slice],
) -> tuple[NDArray[numpy.complex128], tuple[cfdfield_t, ...], int]:
mode_count = len(incoming_modes)
if len(outgoing_modes) != mode_count:
raise ValueError('incoming and outgoing monitor bases must contain the same number of modes')
overlaps = (
*(
compute_overlap_e(
E=mode['E'],
wavenumber=mode['wavenumber'],
dxes=dxes,
axis=axis,
polarity=incoming_polarity,
slices=slices,
)
for mode in incoming_modes
),
*(
compute_overlap_e(
E=mode['E'],
wavenumber=mode['wavenumber'],
dxes=dxes,
axis=axis,
polarity=outgoing_polarity,
slices=slices,
)
for mode in outgoing_modes
),
)
basis_fields = (
*(
expand_e(
E=mode['E'],
wavenumber=mode['wavenumber'],
dxes=dxes,
axis=axis,
polarity=incoming_polarity,
slices=slices,
)
for mode in incoming_modes
),
*(
expand_e(
E=mode['E'],
wavenumber=mode['wavenumber'],
dxes=dxes,
axis=axis,
polarity=outgoing_polarity,
slices=slices,
)
for mode in outgoing_modes
),
)
matrix = numpy.empty((2 * mode_count, 2 * mode_count), dtype=complex)
for row, overlap in enumerate(overlaps):
overlap_vec_conj = vec(overlap).conj()
for col, basis_field in enumerate(basis_fields):
matrix[row, col] = vec(basis_field) @ overlap_vec_conj
return matrix, overlaps, mode_count
def _decompose_port_field(
decomposition: tuple[NDArray[numpy.complex128], tuple[cfdfield_t, ...], int],
field: cfdfield_t,
) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]:
matrix, overlaps, mode_count = decomposition
field_vec = vec(field)
measured = numpy.asarray([field_vec @ vec(overlap).conj() for overlap in overlaps])
coeffs, _resid, _rank, _sing = numpy.linalg.lstsq(matrix, measured, rcond=None)
return coeffs[:mode_count], coeffs[mode_count:]
def _copy_cross_section(
target: NDArray,
cross_section: NDArray,
axis: int,
index: int,
) -> None:
slices: list[slice | int] = [slice(None), slice(None), slice(None), slice(None)]
slices[axis + 1] = index
target[tuple(slices)] = cross_section
def _validate_cross_section(name: str, field: NDArray) -> NDArray:
array = numpy.asarray(field)
if array.ndim != 3 or array.shape[0] != 3:
raise ValueError(f'{name} must have shape (3, n0, n1)')
if not numpy.isfinite(array).all():
raise ValueError(f'{name} must contain only finite values')
return array
def _as_axis_pair(name: str, value: int | Sequence[int]) -> tuple[int, int]:
if isinstance(value, int):
pair = (value, value)
else:
pair = tuple(operator.index(v) for v in value)
if len(pair) != 2:
raise ValueError(f'{name} must be an integer or a two-entry sequence')
if min(pair) < 0:
raise ValueError(f'{name} entries must be non-negative')
return pair
def _pad_cross_section(field: NDArray, padding_cells: tuple[int, int]) -> NDArray:
if padding_cells == (0, 0):
return field
return numpy.pad(
field,
(
(0, 0),
(padding_cells[0], padding_cells[0]),
(padding_cells[1], padding_cells[1]),
),
mode='edge',
)
def _pad_dxes_2d(dxes_2d: dx_lists2_t, padding_cells: tuple[int, int]) -> list[list[NDArray[floating | complexfloating]]]:
padded: list[list[NDArray[floating | complexfloating]]] = [[numpy.empty(0), numpy.empty(0)] for _ in range(2)]
for eh in range(2):
for axis_2d in range(2):
values = numpy.asarray(dxes_2d[eh][axis_2d])
pad = padding_cells[axis_2d]
if pad == 0:
padded[eh][axis_2d] = values.copy()
else:
padded[eh][axis_2d] = numpy.pad(values, (pad, pad), mode='edge')
return padded
def get_s_fdfd(
omega: complex,
dxes: dx_lists_t,
epsilon: fdfield,
*,
left_mode_numbers: Sequence[int],
left_source_slice: Sequence[slice],
left_monitor_slice: Sequence[slice],
right_source_slice: Sequence[slice],
right_monitor_slice: Sequence[slice],
right_mode_numbers: Sequence[int] | None = None,
axis: int = 0,
solve_dxes: dx_lists_t | None = None,
mu: fdfield | None = None,
matrix_solver: Callable[..., Any] | None = None,
matrix_solver_opts: dict[str, Any] | None = None,
source_scales: Sequence[complex] | None = None,
incident_amplitude_atol: float = 1e-12,
return_diagnostics: bool = False,
) -> NDArray[numpy.complex128] | tuple[NDArray[numpy.complex128], dict[str, Any]]:
"""
Extract a two-port guided-mode scattering matrix from local 3D FDFD solves.
The returned S-matrix is ordered as `[[R12, T12], [T21, R21]]`, with
reference planes at the supplied monitor slices. Each column is normalized
by the incident modal amplitude measured on the corresponding input monitor.
"""
if incident_amplitude_atol <= 0:
raise ValueError('incident_amplitude_atol must be positive')
epsilon = numpy.asarray(epsilon)
if epsilon.ndim != 4 or epsilon.shape[0] != 3:
raise ValueError('epsilon must have shape (3, nx, ny, nz)')
if mu is not None:
mu = numpy.asarray(mu)
if mu.shape != epsilon.shape:
raise ValueError('mu must have the same shape as epsilon')
if right_mode_numbers is None:
right_mode_numbers = left_mode_numbers
shape = epsilon.shape[1:]
left_modes_i, right_modes_i, left_source, left_monitor, right_source, right_monitor = _validate_fdfd_port_specs(
left_mode_numbers,
right_mode_numbers,
left_source_slice,
left_monitor_slice,
right_source_slice,
right_monitor_slice,
axis,
shape,
)
n_left = len(left_modes_i)
n_right = len(right_modes_i)
n_total = n_left + n_right
if solve_dxes is None:
solve_dxes = dxes
if source_scales is None:
source_scale_array = numpy.ones((n_total,), dtype=complex)
else:
source_scale_array = numpy.asarray(source_scales, dtype=complex)
if source_scale_array.shape != (n_total,):
raise ValueError('source_scales must have one entry per incident mode')
def solve_modes(mode_numbers: Sequence[int], polarity: int, slices: Sequence[slice]) -> tuple[Waveguide3DMode, ...]:
return tuple(
solve_mode(mode_number, omega=omega, dxes=dxes, axis=axis, polarity=polarity, slices=slices, epsilon=epsilon, mu=mu)
for mode_number in mode_numbers
)
left_source_modes = solve_modes(left_modes_i, +1, left_source)
right_source_modes = solve_modes(right_modes_i, -1, right_source)
left_in_modes = solve_modes(left_modes_i, +1, left_monitor)
left_out_modes = solve_modes(left_modes_i, -1, left_monitor)
right_in_modes = solve_modes(right_modes_i, -1, right_monitor)
right_out_modes = solve_modes(right_modes_i, +1, right_monitor)
left_decomposition = _port_decomposition(
left_in_modes,
left_out_modes,
dxes=dxes,
axis=axis,
incoming_polarity=+1,
outgoing_polarity=-1,
slices=left_monitor,
)
right_decomposition = _port_decomposition(
right_in_modes,
right_out_modes,
dxes=dxes,
axis=axis,
incoming_polarity=-1,
outgoing_polarity=+1,
slices=right_monitor,
)
s = numpy.zeros((n_total, n_total), dtype=complex)
raw_outgoing = numpy.zeros_like(s)
incident_amplitudes = numpy.zeros((n_total,), dtype=complex)
electric_fields: list[cfdfield_t] = []
magnetic_fields: list[cfdfield_t] = []
excitation_modes: list[tuple[int, int, Sequence[slice], Waveguide3DMode]] = [
*((col, +1, left_source, left_source_modes[col]) for col in range(n_left)),
*((n_left + col, -1, right_source, right_source_modes[col]) for col in range(n_right)),
]
for column, polarity, source_slice, source_mode in excitation_modes:
current = compute_source(
E=source_mode['E'],
wavenumber=source_mode['wavenumber'],
omega=omega,
dxes=dxes,
axis=axis,
polarity=polarity,
slices=source_slice,
epsilon=epsilon,
mu=mu,
)
current = cfdfield_t(source_scale_array[column] * current)
solve_kwargs: dict[str, Any] = {
'omega': omega,
'dxes': solve_dxes,
'J': vec(current),
'epsilon': vec(epsilon),
'mu': vec(mu),
'matrix_solver_opts': {} if matrix_solver_opts is None else dict(matrix_solver_opts),
}
if matrix_solver is not None:
solve_kwargs['matrix_solver'] = matrix_solver
e_field = unvec(solvers.generic(**solve_kwargs), shape)
left_in, left_out = _decompose_port_field(left_decomposition, e_field)
right_in, right_out = _decompose_port_field(right_decomposition, e_field)
if column < n_left:
incident = left_in[column]
else:
incident = right_in[column - n_left]
if abs(incident) <= incident_amplitude_atol:
raise ValueError(f'incident modal amplitude for column {column} is zero')
raw_column = numpy.concatenate((left_out, right_out))
incident_amplitudes[column] = incident
raw_outgoing[:, column] = raw_column
s[:, column] = raw_column / incident
if return_diagnostics:
electric_fields.append(cfdfield_t(e_field))
magnetic_fields.append(cfdfield_t(functional.e2h(omega, solve_dxes, mu=mu)(e_field)))
if not return_diagnostics:
return s
guided_powers = numpy.sum(numpy.abs(s) ** 2, axis=0).real
diagnostics: dict[str, Any] = {
'axis': axis,
'left_mode_numbers': left_modes_i,
'right_mode_numbers': right_modes_i,
'left_source_slice': left_source,
'left_monitor_slice': left_monitor,
'right_source_slice': right_source,
'right_monitor_slice': right_monitor,
'left_wavenumbers': numpy.asarray([mode['wavenumber'] for mode in left_out_modes]),
'right_wavenumbers': numpy.asarray([mode['wavenumber'] for mode in right_out_modes]),
'source_scales': source_scale_array,
'electric_fields': tuple(electric_fields),
'magnetic_fields': tuple(magnetic_fields),
'incident_amplitudes': incident_amplitudes,
'raw_outgoing_amplitudes': raw_outgoing,
'guided_powers': guided_powers,
'unresolved_powers': 1.0 - guided_powers,
'left_modes': left_out_modes,
'right_modes': right_out_modes,
}
return s, diagnostics
def get_s_fdfd_from_2d_ports(
left_epsilon: NDArray,
right_epsilon: NDArray,
*,
omega: complex,
dxes_2d: dx_lists2_t,
mode_numbers: Sequence[int],
right_mode_numbers: Sequence[int] | None = None,
axis: int = 0,
left_mu: NDArray | None = None,
right_mu: NDArray | None = None,
prop_dx: float | None = None,
pml_cells: int = 8,
pml_buffer_cells: int = 2,
source_monitor_cells: int = 5,
monitor_interface_cells: int = 6,
transverse_padding_cells: int | Sequence[int] = 0,
transverse_pml_cells: int | Sequence[int] = 0,
pml_epsilon_effective: float = 1.0,
matrix_solver: Callable[..., Any] | None = None,
matrix_solver_opts: dict[str, Any] | None = None,
source_scales: Sequence[complex] | None = None,
return_diagnostics: bool = False,
) -> NDArray[numpy.complex128] | tuple[NDArray[numpy.complex128], dict[str, Any]]:
"""
Extrude two 2D port cross-sections into a short 3D FDFD interface problem.
"""
if axis not in (0, 1, 2):
raise ValueError('axis must be 0, 1, or 2')
if min(pml_cells, pml_buffer_cells, source_monitor_cells, monitor_interface_cells) < 1:
raise ValueError('PML, buffer, source/monitor, and monitor/interface cell counts must be positive')
padding_cells = _as_axis_pair('transverse_padding_cells', transverse_padding_cells)
transverse_pml = _as_axis_pair('transverse_pml_cells', transverse_pml_cells)
left_epsilon = _pad_cross_section(_validate_cross_section('left_epsilon', left_epsilon), padding_cells)
right_epsilon = _pad_cross_section(_validate_cross_section('right_epsilon', right_epsilon), padding_cells)
if left_epsilon.shape != right_epsilon.shape:
raise ValueError('left_epsilon and right_epsilon must have the same shape')
if left_mu is not None:
left_mu = _pad_cross_section(_validate_cross_section('left_mu', left_mu), padding_cells)
if left_mu.shape != left_epsilon.shape:
raise ValueError('left_mu must match left_epsilon')
if right_mu is not None:
right_mu = _pad_cross_section(_validate_cross_section('right_mu', right_mu), padding_cells)
if right_mu.shape != right_epsilon.shape:
raise ValueError('right_mu must match right_epsilon')
if (left_mu is None) != (right_mu is None):
raise ValueError('left_mu and right_mu must be supplied together')
if right_mode_numbers is None:
right_mode_numbers = mode_numbers
dxes_2d_padded = _pad_dxes_2d(dxes_2d, padding_cells)
transverse_shape = left_epsilon.shape[1:]
transverse_axes = tuple(ax for ax in (0, 1, 2) if ax != axis)
if prop_dx is None:
prop_dx = float(numpy.median(numpy.real(numpy.concatenate([numpy.asarray(d) for d in dxes_2d_padded[0]]))))
if prop_dx <= 0:
raise ValueError('prop_dx must be positive')
left_source_i = pml_cells + pml_buffer_cells
left_monitor_i = left_source_i + source_monitor_cells
interface_i = left_monitor_i + monitor_interface_cells
right_monitor_i = interface_i + monitor_interface_cells
right_source_i = right_monitor_i + source_monitor_cells
axis_size = right_source_i + pml_buffer_cells + pml_cells + 1
shape3 = [0, 0, 0]
shape3[axis] = axis_size
shape3[transverse_axes[0]] = transverse_shape[0]
shape3[transverse_axes[1]] = transverse_shape[1]
epsilon = numpy.empty((3, *shape3), dtype=numpy.result_type(left_epsilon, right_epsilon))
mu = None if left_mu is None else numpy.empty((3, *shape3), dtype=numpy.result_type(left_mu, right_mu))
for index in range(axis_size):
cross_epsilon = left_epsilon if index < interface_i else right_epsilon
_copy_cross_section(epsilon, cross_epsilon, axis, index)
if mu is not None:
assert left_mu is not None
assert right_mu is not None
cross_mu = left_mu if index < interface_i else right_mu
_copy_cross_section(mu, cross_mu, axis, index)
dxes: list[list[NDArray[floating | complexfloating]]] = [[numpy.empty(0), numpy.empty(0), numpy.empty(0)] for _ in range(2)]
for eh in range(2):
dxes[eh][axis] = numpy.full(axis_size, prop_dx, dtype=float)
dxes[eh][transverse_axes[0]] = dxes_2d_padded[eh][0].copy()
dxes[eh][transverse_axes[1]] = dxes_2d_padded[eh][1].copy()
solve_dxes: dx_lists_t = [[d.copy() for d in dx] for dx in dxes]
for stretch_axis, thickness in ((axis, pml_cells), (transverse_axes[0], transverse_pml[0]), (transverse_axes[1], transverse_pml[1])):
if thickness == 0:
continue
solve_dxes = scpml.stretch_with_scpml(
solve_dxes,
axis=stretch_axis,
polarity=+1,
omega=float(numpy.real(omega)),
epsilon_effective=pml_epsilon_effective,
thickness=thickness,
)
solve_dxes = scpml.stretch_with_scpml(
solve_dxes,
axis=stretch_axis,
polarity=-1,
omega=float(numpy.real(omega)),
epsilon_effective=pml_epsilon_effective,
thickness=thickness,
)
def port_slice(index: int) -> tuple[slice, slice, slice]:
slices_l = [slice(None), slice(None), slice(None)]
slices_l[axis] = slice(index, index + 1)
return (slices_l[0], slices_l[1], slices_l[2])
result = get_s_fdfd(
omega=omega,
dxes=dxes,
solve_dxes=solve_dxes,
epsilon=epsilon,
mu=mu,
axis=axis,
left_mode_numbers=mode_numbers,
left_source_slice=port_slice(left_source_i),
left_monitor_slice=port_slice(left_monitor_i),
right_mode_numbers=right_mode_numbers,
right_source_slice=port_slice(right_source_i),
right_monitor_slice=port_slice(right_monitor_i),
matrix_solver=matrix_solver,
matrix_solver_opts=matrix_solver_opts,
source_scales=source_scales,
return_diagnostics=return_diagnostics,
)
if not return_diagnostics:
return result
s, diagnostics = result
diagnostics.update({
'interface_index': interface_i,
'pml_cells': pml_cells,
'pml_buffer_cells': pml_buffer_cells,
'source_monitor_cells': source_monitor_cells,
'monitor_interface_cells': monitor_interface_cells,
'transverse_padding_cells': padding_cells,
'transverse_pml_cells': transverse_pml,
'prop_dx': prop_dx,
'epsilon': epsilon,
'mu': mu,
'dxes': dxes,
'solve_dxes': solve_dxes,
})
return s, diagnostics

207
meanas/test/test_eme_numerics.py
View file

@ -5,7 +5,7 @@ import pytest
from scipy import sparse
from ..fdmath import vec
from ..fdfd import eme, waveguide_2d, waveguide_cyl
from ..fdfd import eme, waveguide_2d, waveguide_3d, waveguide_cyl
from ._test_builders import complex_ramp, unit_dxes
from .utils import assert_close
@ -15,7 +15,10 @@ DXES = unit_dxes((2, 2))
WAVENUMBERS_L = numpy.array([1.0, 0.8])
WAVENUMBERS_R = numpy.array([0.9, 0.7])
OMEGA = 1 / 1500
FDFD_OMEGA = 2 * numpy.pi / 16
REAL_DXES = unit_dxes((5, 5))
FDFD_SHAPE_2D = (7, 7)
FDFD_DXES_2D = unit_dxes(FDFD_SHAPE_2D)
def _mode(scale: float) -> tuple[numpy.ndarray, numpy.ndarray]:
@ -481,6 +484,25 @@ def _expected_bragg_reflector_s(n_low: float, n_high: float, periods: int) -> nu
)
def _fdfd_cross_section(core_epsilon: float) -> numpy.ndarray:
epsilon = numpy.ones((3, *FDFD_SHAPE_2D), dtype=float)
epsilon[:, 2:5, 2:5] = core_epsilon
return epsilon
def _small_fdfd_kwargs() -> dict:
return {
'omega': FDFD_OMEGA,
'dxes_2d': FDFD_DXES_2D,
'mode_numbers': [0],
'pml_cells': 4,
'pml_buffer_cells': 2,
'source_monitor_cells': 5,
'monitor_interface_cells': 5,
'matrix_solver_opts': {'atol': 1e-10, 'rtol': 1e-8},
}
def test_get_s_is_near_identity_for_identical_solved_straight_modes() -> None:
mode, wavenumber, _epsilon = _build_straight_mode()
@ -555,6 +577,189 @@ def test_get_s_accepts_complex_leaky_wavenumbers_for_abrupt_interface() -> None:
assert_close(ss, numpy.array([[0.0, 1.0], [1.0, 0.0]], dtype=complex), atol=1e-12, rtol=1e-12)
def test_get_s_fdfd_from_2d_ports_identical_guides_are_nearly_transmissive() -> None:
epsilon = _fdfd_cross_section(2.25)
ss = waveguide_3d.get_s_fdfd_from_2d_ports(epsilon, epsilon, **_small_fdfd_kwargs())
assert ss.shape == (2, 2)
assert numpy.isfinite(ss).all()
assert abs(ss[0, 0]) < 0.03
assert abs(ss[1, 1]) < 0.04
assert abs(abs(ss[1, 0]) - 1.0) < 0.02
assert abs(abs(ss[0, 1]) - 1.0) < 0.03
def test_get_s_fdfd_from_2d_ports_mismatch_returns_finite_diagnostics() -> None:
left = _fdfd_cross_section(2.25)
right = _fdfd_cross_section(1.5)
ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(left, right, return_diagnostics=True, **_small_fdfd_kwargs())
assert ss.shape == (2, 2)
assert numpy.isfinite(ss).all()
assert len(diagnostics['electric_fields']) == 2
assert len(diagnostics['magnetic_fields']) == 2
assert diagnostics['electric_fields'][0].shape == (3, 33, *FDFD_SHAPE_2D)
assert numpy.isfinite(diagnostics['electric_fields'][0]).all()
assert numpy.isfinite(diagnostics['magnetic_fields'][0]).all()
assert numpy.isfinite(diagnostics['guided_powers']).all()
assert numpy.isfinite(diagnostics['unresolved_powers']).all()
assert numpy.isfinite(diagnostics['raw_outgoing_amplitudes']).all()
assert numpy.isfinite(diagnostics['incident_amplitudes']).all()
assert (numpy.abs(diagnostics['incident_amplitudes']) > 1e-3).all()
assert_close(diagnostics['unresolved_powers'], 1.0 - diagnostics['guided_powers'])
assert diagnostics['guided_powers'][0] < 1.05
def test_get_s_fdfd_normalization_is_source_scale_invariant() -> None:
left = _fdfd_cross_section(2.25)
right = _fdfd_cross_section(1.5)
kwargs = _small_fdfd_kwargs()
base, base_diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(left, right, return_diagnostics=True, **kwargs)
scaled, scaled_diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(left, right, source_scales=[2.0, 0.3], return_diagnostics=True, **kwargs)
assert_close(numpy.abs(scaled), numpy.abs(base), atol=1e-6, rtol=1e-6)
assert abs(abs(scaled_diagnostics['incident_amplitudes'][0] / base_diagnostics['incident_amplitudes'][0]) - 2.0) < 1e-6
assert abs(abs(scaled_diagnostics['incident_amplitudes'][1] / base_diagnostics['incident_amplitudes'][1]) - 0.3) < 1e-6
def test_get_s_fdfd_from_2d_ports_diagnostics_order_generated_slices() -> None:
epsilon = _fdfd_cross_section(2.25)
_ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(epsilon, epsilon, return_diagnostics=True, **_small_fdfd_kwargs())
left_source = diagnostics['left_source_slice'][0].start
left_monitor = diagnostics['left_monitor_slice'][0].start
interface_index = diagnostics['interface_index']
right_monitor = diagnostics['right_monitor_slice'][0].start
right_source = diagnostics['right_source_slice'][0].start
assert left_source == 6
assert left_monitor == 11
assert interface_index == 16
assert right_monitor == 21
assert right_source == 26
assert left_source < left_monitor < interface_index < right_monitor < right_source
def test_get_s_fdfd_from_2d_ports_transverse_padding_and_pml_updates_diagnostics() -> None:
epsilon = _fdfd_cross_section(2.25)
ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(
epsilon,
epsilon,
transverse_padding_cells=(2, 3),
transverse_pml_cells=(1, 2),
return_diagnostics=True,
**_small_fdfd_kwargs(),
)
assert ss.shape == (2, 2)
assert diagnostics['epsilon'].shape == (3, 33, 11, 13)
assert diagnostics['transverse_padding_cells'] == (2, 3)
assert diagnostics['transverse_pml_cells'] == (1, 2)
assert numpy.iscomplexobj(diagnostics['solve_dxes'][0][1])
assert numpy.iscomplexobj(diagnostics['solve_dxes'][0][2])
assert numpy.abs(numpy.imag(diagnostics['solve_dxes'][0][1][0])) > 0
assert numpy.abs(numpy.imag(diagnostics['solve_dxes'][0][2][0])) > 0
def test_get_s_fdfd_from_2d_ports_supports_nondefault_axis() -> None:
epsilon = _fdfd_cross_section(2.25)
ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(epsilon, epsilon, axis=1, return_diagnostics=True, **_small_fdfd_kwargs())
assert ss.shape == (2, 2)
assert diagnostics['epsilon'].shape == (3, FDFD_SHAPE_2D[0], 33, FDFD_SHAPE_2D[1])
assert diagnostics['left_source_slice'][1].start == 6
assert numpy.isfinite(ss).all()
def test_get_s_fdfd_rejects_bad_port_slice_order() -> None:
epsilon = numpy.ones((3, 12, 5, 5), dtype=float)
dxes = [[numpy.ones(12), numpy.ones(5), numpy.ones(5)] for _ in range(2)]
with pytest.raises(ValueError, match='ordered'):
waveguide_3d.get_s_fdfd(
FDFD_OMEGA,
dxes,
epsilon,
left_mode_numbers=[0],
left_source_slice=(slice(5, 6), slice(None), slice(None)),
left_monitor_slice=(slice(4, 5), slice(None), slice(None)),
right_source_slice=(slice(10, 11), slice(None), slice(None)),
right_monitor_slice=(slice(8, 9), slice(None), slice(None)),
)
def test_get_s_fdfd_rejects_bad_source_scale_count_before_solving() -> None:
epsilon = numpy.ones((3, 12, 5, 5), dtype=float)
dxes = [[numpy.ones(12), numpy.ones(5), numpy.ones(5)] for _ in range(2)]
with pytest.raises(ValueError, match='source_scales'):
waveguide_3d.get_s_fdfd(
FDFD_OMEGA,
dxes,
epsilon,
left_mode_numbers=[0],
left_source_slice=(slice(2, 3), slice(None), slice(None)),
left_monitor_slice=(slice(4, 5), slice(None), slice(None)),
right_mode_numbers=[0, 1],
right_source_slice=(slice(10, 11), slice(None), slice(None)),
right_monitor_slice=(slice(8, 9), slice(None), slice(None)),
source_scales=[1.0, 1.0],
)
def test_get_s_fdfd_from_2d_ports_rejects_invalid_generated_geometry() -> None:
epsilon = _fdfd_cross_section(2.25)
with pytest.raises(ValueError, match='must be positive'):
waveguide_3d.get_s_fdfd_from_2d_ports(epsilon, epsilon, **{**_small_fdfd_kwargs(), 'pml_cells': 0})
def test_get_s_fdfd_rejects_noninteger_mode_numbers_before_solving() -> None:
epsilon = numpy.ones((3, 12, 5, 5), dtype=float)
dxes = [[numpy.ones(12), numpy.ones(5), numpy.ones(5)] for _ in range(2)]
with pytest.raises(ValueError, match='mode_numbers'):
waveguide_3d.get_s_fdfd(
FDFD_OMEGA,
dxes,
epsilon,
left_mode_numbers=[0.5], # type: ignore[list-item]
left_source_slice=(slice(2, 3), slice(None), slice(None)),
left_monitor_slice=(slice(4, 5), slice(None), slice(None)),
right_source_slice=(slice(10, 11), slice(None), slice(None)),
right_monitor_slice=(slice(8, 9), slice(None), slice(None)),
)
def test_get_s_fdfd_from_2d_ports_preserves_complex_generated_inputs(monkeypatch: pytest.MonkeyPatch) -> None:
left = _fdfd_cross_section(2.25)
right = _fdfd_cross_section(1.5).astype(complex)
right += 0.1j
dxes_2d = [
[FDFD_DXES_2D[0][0].astype(complex) + 0.2j, FDFD_DXES_2D[0][1].astype(complex) + 0.3j],
[FDFD_DXES_2D[1][0].astype(complex) + 0.4j, FDFD_DXES_2D[1][1].astype(complex) + 0.5j],
]
def fake_get_s_fdfd(*args, **kwargs) -> tuple[numpy.ndarray, dict]:
return numpy.eye(2, dtype=complex), {}
monkeypatch.setattr(waveguide_3d, 'get_s_fdfd', fake_get_s_fdfd)
kwargs = {**_small_fdfd_kwargs(), 'dxes_2d': dxes_2d}
_ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(left, right, return_diagnostics=True, **kwargs)
assert numpy.iscomplexobj(diagnostics['epsilon'])
assert diagnostics['epsilon'][0, diagnostics['interface_index'], 0, 0].imag == pytest.approx(0.1)
assert numpy.iscomplexobj(diagnostics['dxes'][0][1])
assert diagnostics['dxes'][0][1][0].imag == pytest.approx(0.2)
def test_quarter_wave_matching_layer_is_nearly_reflectionless_at_design_frequency() -> None:
"""
A single quarter-wave matching layer with