From 37ed9d10bcf3b9cbae4ad0a43c0040d04e41ba55 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Tue, 14 Jul 2026 00:21:12 -0700 Subject: [PATCH] [fdfd.waveguide_3d] add FDFD-based S-matrix solve --- meanas/fdfd/waveguide_3d.py | 560 ++++++++++++++++++++++++++++++- meanas/test/test_eme_numerics.py | 207 +++++++++++- 2 files changed, 762 insertions(+), 5 deletions(-) diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py index c77c8d4..3f7a29b 100644 --- a/meanas/fdfd/waveguide_3d.py +++ b/meanas/fdfd/waveguide_3d.py @@ -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 diff --git a/meanas/test/test_eme_numerics.py b/meanas/test/test_eme_numerics.py index 0d28692..45ad7cc 100644 --- a/meanas/test/test_eme_numerics.py +++ b/meanas/test/test_eme_numerics.py @@ -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