diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py index 958d1d0..a0949e3 100644 --- a/meanas/fdfd/waveguide_2d.py +++ b/meanas/fdfd/waveguide_2d.py @@ -801,77 +801,40 @@ def sensitivity( wavenumber: complex, omega: complex, dxes: dx_lists2_t, - epsilon: vfdslice, - mu: vfdslice | None = None, + epsilon: vfdslice | vcfdslice, + mu: vfdslice | vcfdslice | None = None, ) -> vcfdslice_t: r""" - Given a waveguide structure (`dxes`, `epsilon`, `mu`) and mode fields - (`e_norm`, `h_norm`, `wavenumber`, `omega`), calculates the sensitivity of the wavenumber - $\beta$ to changes in the dielectric structure $\epsilon$. + Calculate the analytic sensitivity of an isolated mode's beta to epsilon. - The output is a vector of the same size as `vec(epsilon)`, with each element specifying the - sensitivity of `wavenumber` to changes in the corresponding element in `vec(epsilon)`, i.e. + For the reciprocal diagonal waveguide problem ``A_E E_xy = beta**2 E_xy``, + the Yee-weighted Lorentz left eigenvector is the unconjugated row - $$sens_{i} = \frac{\partial\beta}{\partial\epsilon_i}$$ + $$ L^T = [H_y \Delta x_e \Delta y_h,\; -H_x \Delta x_h \Delta y_e]. $$ - An adjoint approach is used to calculate the sensitivity; the derivation is provided here: + Differentiating the eigenproblem and projecting with this left vector gives - Starting with the eigenvalue equation - - $$\beta^2 E_{xy} = A_E E_{xy}$$ - - where $A_E$ is the waveguide operator from `operator_e()`, and $E_{xy} = \begin{bmatrix} E_x \\ - E_y \end{bmatrix}$, - we can differentiate with respect to one of the $\epsilon$ elements (i.e. at one Yee grid point), $\epsilon_i$: - - $$ - (2 \beta) \partial_{\epsilon_i}(\beta) E_{xy} + \beta^2 \partial_{\epsilon_i} E_{xy} - = \partial_{\epsilon_i}(A_E) E_{xy} + A_E \partial_{\epsilon_i} E_{xy} - $$ - - We then multiply by $H_{yx}^\star = \begin{bmatrix}H_y^\star \\ -H_x^\star \end{bmatrix}$ from the left: - - $$ - (2 \beta) \partial_{\epsilon_i}(\beta) H_{yx}^\star E_{xy} + \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy} - = H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} + H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} - $$ - - On a graded Yee grid, the left projection includes the quadrature weights: - $H_y^\star \Delta x_e \Delta y_h$ and $-H_x^\star \Delta x_h \Delta y_e$. - These weights multiply the left vector, before it acts on any derivative operators. - With these weights understood, $H_{yx}^\star$ is a left-eigenvector of $A_E$. This can be verified by inspecting - the form of `operator_h` ($A_H$) and comparing its conjugate transpose to `operator_e` ($A_E$). Also, note - $H_{yx}^\star \cdot E_{xy} = H^\star \times E$ recalls the mode orthogonality relation. See doi:10.5194/ars-9-85-201 - for a similar approach. Therefore, - - $$ - H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} = \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy} - $$ - - and we can simplify to - - $$ - \partial_{\epsilon_i}(\beta) - = \frac{1}{2 \beta} \frac{H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} }{H_{yx}^\star E_{xy}} - $$ - - This expression can be quickly calculated for all $i$ by writing out the various terms of - $\partial_{\epsilon_i} A_E$ and recognizing that the vector-matrix-vector products (i.e. scalars) - $sens_i = \vec{v}_{left} \partial_{\epsilon_i} (\epsilon_{xyz}) \vec{v}_{right}$, indexed by $i$, can be expressed as - elementwise multiplications $\vec{sens} = \vec{v}_{left} \star \vec{v}_{right}$ + $$ \frac{\partial\beta}{\partial\epsilon_i} + = \frac{L^T (\partial_{\epsilon_i} A_E) E_{xy}} + {2\beta L^T E_{xy}}. $$ + This pairing is bilinear: neither H nor the coordinate weights are + conjugated. It therefore applies to complex PML stretches and reciprocal + lossy materials as well as lossless grids. Quadrature weights multiply the + left projection before derivative operators act. Overall field scaling + cancels between numerator and denominator. Args: - e_norm: Normalized, vectorized E_xyz field for the mode. E.g. as returned by `normalized_fields_e`. - h_norm: Normalized, vectorized H_xyz field for the mode. E.g. as returned by `normalized_fields_e`. - wavenumber: Propagation constant for the mode. The z-axis is assumed to be continuous (i.e. without numerical dispersion). - omega: The angular frequency of the system. - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) - epsilon: Vectorized dielectric constant grid - mu: Vectorized magnetic permeability grid (default 1 everywhere) + e_norm: Vectorized E_xyz eigenfield for the mode. + h_norm: Matching vectorized H_xyz eigenfield at the same phase reference. + wavenumber: Nonzero beta of the continuous-z eigenproblem. + omega: Angular frequency of the system. + dxes: Two-dimensional Yee widths, including any complex PML stretches. + epsilon: Vectorized diagonal permittivity. + mu: Optional vectorized diagonal permeability; defaults to one. Returns: - Sparse matrix representation of the operator. + Vector ``d beta / d epsilon_i`` in the same component order as epsilon. """ if mu is None: mu = numpy.ones_like(epsilon) @@ -891,12 +854,12 @@ def sensitivity( ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2]) hx, hy, _hz = numpy.split(h_norm, 3) # The Yee-weighted H projection is the left eigenvector on a graded grid. - hv_yx_conj = numpy.conj(numpy.concatenate([hy, -hx])) * numpy.concatenate([da_exxhyy, da_eyyhxx]) + hv_yx = numpy.concatenate([hy, -hx]) * numpy.concatenate([da_exxhyy, da_eyyhxx]) - sens_xy1 = (hv_yx_conj @ (omega * omega * mu_yx)) * ev_xy - sens_xy2 = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby))) * ev_xy - sens_z = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ (-eps_z_inv * eps_z_inv)) * (sparse.hstack((Dbx, Dby)) @ eps_xy @ ev_xy) - norm = hv_yx_conj @ ev_xy + sens_xy1 = (hv_yx @ (omega * omega * mu_yx)) * ev_xy + sens_xy2 = (hv_yx @ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby))) * ev_xy + sens_z = (hv_yx @ sparse.vstack((Dfx, Dfy)) @ (-eps_z_inv * eps_z_inv)) * (sparse.hstack((Dbx, Dby)) @ eps_xy @ ev_xy) + norm = hv_yx @ ev_xy sens_tot = numpy.concatenate([sens_xy1 + sens_xy2, sens_z]) / (2 * wavenumber * norm) return vcfdslice_t(sens_tot) diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py index 8955a1a..15efebc 100644 --- a/meanas/fdfd/waveguide_3d.py +++ b/meanas/fdfd/waveguide_3d.py @@ -143,6 +143,9 @@ def mode_from_2d( ``polarity`` is +1 toward increasing indices and -1 toward decreasing indices. ``normalization`` follows waveguide_2d.normalized_fields_e. + Lossless modes at or beyond the longitudinal grid band edge cannot be + power-normalized. Lorentz normalization can represent these modes + beyond that edge and selects the decaying branch for real beta. """ if polarity not in (-1, 1): raise ValueError('polarity must be -1 or +1') @@ -154,7 +157,12 @@ def mode_from_2d( raise ValueError('the 2D eigenvector and wavenumber must be finite and match the cross-section') wavenumber_2d = complex(wavenumber_2d) # Correct wavenumber to account for numerical dispersion. - wavenumber = 2 / dx_prop * numpy.arcsin(wavenumber_2d * dx_prop / 2) + argument = wavenumber_2d * dx_prop / 2 + if numpy.isreal(argument) and abs(argument.real) >= 1: + if normalization == 'power': + raise ValueError('mode is at or beyond the longitudinal grid band edge; refine the propagation grid') + argument = complex(argument.real, -0.0) + wavenumber = 2 / dx_prop * numpy.arcsin(argument) shape = [d.size for d in args_2d['dxes'][0]] ve, vh = waveguide_2d.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, prop_phase=dx_prop * wavenumber, normalization=normalization, **args_2d) @@ -372,12 +380,18 @@ def expand_e( $$ to each copied slice. + + Transverse electric samples use forward intervals dx_e[i] between + indices i and i+1. The longitudinal electric component uses backward + intervals dx_h[i] between indices i-1 and i on its staggered grid. """ slices = tuple(slices) # Determine phase factors for parallel slices a_shape = numpy.roll([1, -1, 1, 1], axis) - a_E = numpy.real(dxes[0][axis]).cumsum() + # dx_e[i] spans transverse electric samples i -> i+1. The longitudinal + # component lies on the complementary grid, with backward intervals dx_h. + a_E = numpy.concatenate(([0.0], numpy.real(dxes[0][axis][:-1]).cumsum())) r_E = a_E - a_E[slices[axis]] iphi = polarity * -1j * wavenumber phase_E = numpy.exp(iphi * r_E).reshape(a_shape) @@ -392,6 +406,9 @@ def expand_e( slices_in = (slice(None), *slices) E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in] + a_long = numpy.real(dxes[1][axis]).cumsum() + phase_long = numpy.exp(iphi * (a_long - a_long[slices[axis]])).reshape(a_shape[1:]) + E_expanded[(axis, *slices_exp_l)] = phase_long * E[(axis, *slices)] return cfdfield_t(E_expanded) @@ -550,10 +567,41 @@ def _decompose_port_field( 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) + coeffs, _resid, rank, _sing = numpy.linalg.lstsq(matrix, measured, rcond=None) + if rank != matrix.shape[1]: + raise numpy.linalg.LinAlgError('monitor modal basis is rank deficient') return coeffs[:mode_count], coeffs[mode_count:] +def _port_power_gram( + modes: Sequence[Waveguide3DMode], + dxes: dx_lists_t, + axis: int, + polarity: int, + slices: Sequence[slice], + ) -> NDArray[numpy.complex128]: + """Hermitian flux metric, including interference between nonorthogonal modes.""" + u, v = (axis + 1) % 3, (axis + 2) % 3 + + def width(family: int, component: int) -> NDArray: + values = numpy.real(dxes[family][component][slices[component]]) + shape = [1, 1, 1] + shape[component] = values.size + return values.reshape(shape) + + area_uv = width(0, u) * width(1, v) + area_vu = width(1, u) * width(0, v) + dx = float(numpy.real(dxes[0][axis][slices[axis]]).item()) + gram = numpy.empty((len(modes), len(modes)), dtype=complex) + for row, test in enumerate(modes): + h = test['H'][(slice(None), *slices)].conj() + for column, mode in enumerate(modes): + e = mode['E'][(slice(None), *slices)] + phase = numpy.exp(-1j * polarity * mode['wavenumber'] * dx) + gram[row, column] = 0.5 * polarity * phase * numpy.sum(e[u] * h[v] * area_uv - e[v] * h[u] * area_vu) + return 0.5 * (gram + gram.conj().T) + + def _copy_cross_section( target: NDArray, cross_section: NDArray, @@ -637,8 +685,11 @@ def get_s_fdfd( Extract a two-port guided-mode scattering matrix from local 3D FDFD solves. The returned S-matrix is ordered as `[[R_L, T_R_to_L], [T_L_to_R, R_R]]`, with - reference planes at the supplied monitor slices. Each column is normalized - by the incident modal amplitude measured on the corresponding input monitor. + reference planes at the supplied monitor slices. Calibration solves + ``S @ incoming_amplitudes = raw_outgoing_amplitudes`` using every measured + incoming coefficient on both ports, including mixed or degenerate modes. + Diagnostic guided powers use the full input/output flux Gram matrices, + retaining interference between nonorthogonal modal basis vectors. """ if incident_amplitude_atol <= 0: raise ValueError('incident_amplitude_atol must be positive') @@ -710,8 +761,8 @@ def get_s_fdfd( slices=right_monitor, ) - s = numpy.zeros((n_total, n_total), dtype=complex) - raw_outgoing = numpy.zeros_like(s) + raw_outgoing = numpy.zeros((n_total, n_total), dtype=complex) + incoming_amplitudes = numpy.zeros_like(raw_outgoing) incident_amplitudes = numpy.zeros((n_total,), dtype=complex) electric_fields: list[cfdfield_t] = [] magnetic_fields: list[cfdfield_t] = [] @@ -747,26 +798,39 @@ def get_s_fdfd( 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: + incoming_column = numpy.concatenate((left_in, right_in)) + if numpy.linalg.norm(incoming_column) <= 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 + incident_amplitudes[column] = incoming_column[column] + incoming_amplitudes[:, column] = incoming_column 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 numpy.isfinite(incoming_amplitudes).all() or not numpy.isfinite(raw_outgoing).all(): + raise ValueError('modal calibration requires finite incoming and outgoing amplitudes') + try: + s = numpy.linalg.solve(incoming_amplitudes.T, raw_outgoing.T).T + except numpy.linalg.LinAlgError as error: + raise numpy.linalg.LinAlgError('incoming excitation matrix is singular') from error + if not return_diagnostics: return s - guided_powers = numpy.sum(numpy.abs(s) ** 2, axis=0).real + zero = numpy.zeros((n_left, n_right), dtype=complex) + input_power_gram = numpy.block([ + [_port_power_gram(left_in_modes, dxes, axis, +1, left_monitor), zero], + [zero.T, _port_power_gram(right_in_modes, dxes, axis, -1, right_monitor)], + ]) + output_power_gram = numpy.block([ + [_port_power_gram(left_out_modes, dxes, axis, -1, left_monitor), zero], + [zero.T, _port_power_gram(right_out_modes, dxes, axis, +1, right_monitor)], + ]) + guided_powers = numpy.diag(s.conj().T @ output_power_gram @ s).real / numpy.diag(input_power_gram).real diagnostics: dict[str, Any] = { 'axis': axis, 'left_mode_numbers': left_modes_i, @@ -781,6 +845,9 @@ def get_s_fdfd( 'electric_fields': tuple(electric_fields), 'magnetic_fields': tuple(magnetic_fields), 'incident_amplitudes': incident_amplitudes, + 'incoming_amplitudes': incoming_amplitudes, + 'incoming_power_gram': input_power_gram, + 'outgoing_power_gram': output_power_gram, 'raw_outgoing_amplitudes': raw_outgoing, 'guided_powers': guided_powers, 'unresolved_powers': 1.0 - guided_powers, diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index 0d1d4d7..91a109b 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -86,7 +86,7 @@ large-radius-compatible identities $$ E_z = -\frac{1}{\imath \beta} T_a \epsilon_{zz}^{-1} +\frac{1}{\imath \beta} \epsilon_{zz}^{-1} \begin{bmatrix} \hat{\partial}_r T_b \epsilon_{rr} & \hat{\partial}_y T_a \epsilon_{yy} @@ -304,6 +304,11 @@ def linear_wavenumbers( Calculate linear wavenumbers (1/distance) based on angular wavenumbers (1/rad) and the mode's energy distribution. + The effective radius is weighted by transverse electric energy, using the + physical Yee cell areas for each component. Radial E uses dx_e * dy_h; + transverse E uses dx_h * dy_e. Real coordinate widths define these areas + and radii when complex PML stretching is present. + Args: e_xys: Vectorized mode fields with shape (num_modes, 2 * x *y) angular_wavenumbers: Wavenumbers assuming fields have theta-dependence of @@ -321,14 +326,19 @@ def linear_wavenumbers( shape2d = (len(dxes[0][0]), len(dxes[0][1])) epsilon2d = unvec(epsilon, shape2d)[:2] - ra = rmin + numpy.cumsum(dxes[0][0]) - rb = rmin + dxes[0][0] / 2.0 + numpy.concatenate(( - numpy.zeros(1, dtype=dxes[1][0].dtype), - numpy.cumsum(dxes[1][0][:-1]), + physical_widths = [[numpy.real(dx) for dx in family] for family in dxes] + ra = rmin + numpy.cumsum(physical_widths[0][0]) + rb = rmin + physical_widths[0][0] / 2.0 + numpy.concatenate(( + numpy.zeros(1), + numpy.cumsum(physical_widths[1][0][:-1]), )) + areas = numpy.stack(( + numpy.outer(physical_widths[0][0], physical_widths[1][1]), + numpy.outer(physical_widths[1][0], physical_widths[0][1]), + )) for ii in range(angular_wavenumbers.size): efield = unvec(e_xys[ii], shape2d, 2) - energy = numpy.real((efield * efield.conj()) * epsilon2d) + energy = numpy.real((efield * efield.conj()) * epsilon2d) * areas er_energy_vs_r = energy[0].sum(axis=1) ey_energy_vs_r = energy[1].sum(axis=1) energy_vs_r = er_energy_vs_r + ey_energy_vs_r @@ -405,7 +415,7 @@ def exy2e( n_pts = dxes[0][0].size * dxes[0][1].size exy2ez = ( - Ta @ epsilon_z_inv + epsilon_z_inv @ sparse.hstack((Dbx @ Tb @ epsilon_x, Dby @ Ta @ epsilon_y)) / (1j * wavenumber) diff --git a/meanas/test/test_eme_numerics.py b/meanas/test/test_eme_numerics.py index b2c078d..7232e73 100644 --- a/meanas/test/test_eme_numerics.py +++ b/meanas/test/test_eme_numerics.py @@ -2,6 +2,7 @@ from typing import cast import numpy import pytest +from scipy.sparse.linalg import spsolve from ..fdmath import vec from ..fdfd import eme, waveguide_2d, waveguide_3d, waveguide_cyl @@ -547,6 +548,50 @@ def test_get_s_fdfd_from_2d_ports_identical_guides_are_nearly_transmissive() -> assert abs(abs(ss[0, 1]) - 1.0) < 0.03 +def test_multimode_fdfd_calibrates_the_full_incoming_matrix(monkeypatch) -> None: + # Exactly degenerate polarizations have different source/monitor bases. + # A uniform guide must still have zero interface reflection and conserve flux. + original_solve = waveguide_2d.solve_mode + seed = 0 + + def seeded_solve(*args, **kwargs): + nonlocal seed + seed += 1 + return original_solve(*args, **kwargs, rng=seed) + + decompositions = [] + original_decomposition = waveguide_3d._port_decomposition + + def capture_decomposition(*args, **kwargs): + result = original_decomposition(*args, **kwargs) + decompositions.append(result) + return result + + monkeypatch.setattr(waveguide_2d, 'solve_mode', seeded_solve) + monkeypatch.setattr(waveguide_3d, '_port_decomposition', capture_decomposition) + epsilon = numpy.ones((3, 2, 2)) + widths = [[numpy.full(2, 0.5), numpy.full(2, 0.5)] for _ in range(2)] + ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports( + epsilon, epsilon, omega=0.5, dxes_2d=widths, mode_numbers=[0, 1], + pml_cells=4, pml_buffer_cells=1, source_monitor_cells=2, monitor_interface_cells=2, + source_scales=[2, 0.3, 1j, 0.7j], matrix_solver=spsolve, return_diagnostics=True, + ) + incoming, outgoing = [], [] + for field in diagnostics['electric_fields']: + li, lo = waveguide_3d._decompose_port_field(decompositions[0], field) + ri, ro = waveguide_3d._decompose_port_field(decompositions[1], field) + incoming.append(numpy.concatenate((li, ri))) + outgoing.append(numpy.concatenate((lo, ro))) + incoming, outgoing = numpy.column_stack(incoming), numpy.column_stack(outgoing) + assert numpy.linalg.norm(incoming - numpy.diag(numpy.diag(incoming))) > 0.1 + assert_close(ss @ incoming, outgoing, atol=1e-11, rtol=1e-11) + assert_close(diagnostics['incoming_amplitudes'], incoming, atol=1e-12) + assert_close(ss[:2, :2], 0, atol=1e-11) + assert_close(ss[2:, 2:], 0, atol=1e-11) + assert_close(ss.conj().T @ diagnostics['outgoing_power_gram'] @ ss, diagnostics['incoming_power_gram'], atol=1e-11) + assert_close(diagnostics['guided_powers'], 1, atol=1e-11) + + 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) diff --git a/meanas/test/test_waveguide_2d_modes.py b/meanas/test/test_waveguide_2d_modes.py index 3edbae8..f8af879 100644 --- a/meanas/test/test_waveguide_2d_modes.py +++ b/meanas/test/test_waveguide_2d_modes.py @@ -1261,6 +1261,7 @@ def test_richer_complex_basis_reduces_driven_fdfd_monitor_residual() -> None: dxes_2d=dxes_2d, mode_numbers=[0], axis=2, + prop_dx=0.15, # Resolve both port modes below the longitudinal band edge. pml_cells=4, pml_buffer_cells=2, source_monitor_cells=4, diff --git a/meanas/test/test_waveguide_2d_numerics.py b/meanas/test/test_waveguide_2d_numerics.py index d590be6..5a244c5 100644 --- a/meanas/test/test_waveguide_2d_numerics.py +++ b/meanas/test/test_waveguide_2d_numerics.py @@ -1,10 +1,11 @@ import numpy +import pytest from numpy.linalg import norm from numpy.testing import assert_allclose from scipy import linalg from ..fdmath import vec -from ..fdfd import waveguide_2d +from ..fdfd import waveguide_2d, scpml OMEGA = 1 / 1500 @@ -121,7 +122,8 @@ def test_waveguide_2d_sensitivity_matches_finite_difference() -> None: assert 0.4 < ratio < 1.8 -def test_sensitivity_matches_central_differences_on_graded_yee_mesh() -> None: +@pytest.mark.parametrize('absorption', ['none', 'pml', 'material']) +def test_sensitivity_matches_central_differences_on_graded_yee_mesh(absorption) -> None: rng = numpy.random.default_rng(4) shape = (4, 3) omega = 0.7 @@ -129,6 +131,10 @@ def test_sensitivity_matches_central_differences_on_graded_yee_mesh() -> None: magnetic_widths = [rng.uniform(0.2, 1.2, n) for n in shape] electric_widths = [0.5 * (dx + numpy.roll(dx, -1)) for dx in magnetic_widths] dxes = [electric_widths, magnetic_widths] + if absorption == 'pml': + dxes = scpml.stretch_with_scpml(dxes, axis=0, polarity=-1, omega=omega, thickness=1, max_imag_ratio=0.3) + elif absorption == 'material': + epsilon = epsilon * (1 - 0.15j * rng.uniform(0.5, 1.5, epsilon.size)) values, vectors = linalg.eig(waveguide_2d.operator_e(omega, dxes, epsilon).toarray()) selected = numpy.argmax(values.real) beta = numpy.sqrt(values[selected]) diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py index e5a45a6..3c14d2f 100644 --- a/meanas/test/test_waveguide_mode_helpers.py +++ b/meanas/test/test_waveguide_mode_helpers.py @@ -153,8 +153,8 @@ def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent( line_selection[axis] = slice(None) values = expanded[component][tuple(line_selection)] ratios = values[1:] / values[:-1] - phase_positions = numpy.real(dxes[0][axis]).cumsum() - expected_ratios = numpy.exp(-1j * result['wavenumber'] * numpy.diff(phase_positions)) + intervals = numpy.real(dxes[1][axis][1:] if component == axis else dxes[0][axis][:-1]) + expected_ratios = numpy.exp(-1j * result['wavenumber'] * intervals) numpy.testing.assert_allclose(ratios, expected_ratios, rtol=1e-6, atol=1e-9) @@ -422,7 +422,7 @@ def test_waveguide_cyl_reconstruction_is_finite_radius_consistent(use_mu: bool) eps_r, eps_y, eps_z = numpy.split(epsilon, 3) mu_r, mu_y, mu_z = numpy.split(mu_full, 3) - ez_rhs = ta * ( + ez_rhs = ( (Dbx @ (tb * eps_r * er)) + (Dby @ (ta * eps_y * ey)) ) / eps_z @@ -469,8 +469,79 @@ def test_waveguide_cyl_linear_wavenumbers_use_component_radii() -> None: rmin=rmin, ) - numpy.testing.assert_allclose(er_linear[0], angular_wavenumber[0] / rb.mean()) - numpy.testing.assert_allclose(ey_linear[0], angular_wavenumber[0] / ra.mean()) + numpy.testing.assert_allclose(er_linear[0], angular_wavenumber[0] / numpy.average(rb, weights=dxes[0][0])) + numpy.testing.assert_allclose(ey_linear[0], angular_wavenumber[0] / numpy.average(ra, weights=dxes[1][0])) + + +@pytest.mark.parametrize('use_mu', [False, True]) +def test_cylindrical_modes_obey_the_full_ampere_curl(use_mu) -> None: + shape, omega, rmin = (4, 3), 0.7, 2.0 + dxes = [[numpy.full(n, 0.4) for n in shape] for _ in range(2)] + rng = numpy.random.default_rng(4) + epsilon = vec(1 + rng.uniform(0, 3, (3, *shape))) + mu = vec(1 + rng.uniform(0, 2, (3, *shape))) if use_mu else numpy.ones_like(epsilon) + matrix = waveguide_cyl.cylindrical_operator(omega, dxes, epsilon, rmin, mu).toarray() + values, vectors = linalg.eig(matrix) + selected = numpy.argmax(values.real) + beta = numpy.sqrt(values[selected]) + electric = waveguide_cyl.exy2e(beta * rmin, dxes, rmin, epsilon) @ vectors[:, selected] + magnetic = waveguide_cyl.e2h(beta * rmin, omega, dxes, rmin, mu) @ electric + ta, tb = (metric.diagonal() for metric in waveguide_cyl.dxes2T(dxes, rmin)) + dbx, dby = deriv_back(dxes[1]) + hr, hy, hz = numpy.split(magnetic, 3) + # Physical cylindrical Ampere curl, independently of electric reconstruction. + curl_h = numpy.concatenate(( + 1j * beta * hy / tb + dby @ hz, + -1j * beta * hr / ta - (dbx @ (tb * hz)) / ta, + -dby @ hr + dbx @ hy, + )) + numpy.testing.assert_allclose(curl_h, 1j * omega * epsilon * electric, atol=1e-11, rtol=1e-10) + + +def test_cylindrical_mean_radius_is_invariant_under_constant_density_refinement() -> None: + for widths in ([1.0, 3.0], [1.0, 1.0, 1.0, 1.0]): + n = len(widths) + dxes = [[numpy.array(widths), numpy.ones(1)] for _ in range(2)] + electric_xy = numpy.concatenate((numpy.ones(n), numpy.zeros(n)))[None, :].astype(complex) + beta = waveguide_cyl.linear_wavenumbers(electric_xy, [12.0], numpy.ones(3 * n), dxes, rmin=10.0) + numpy.testing.assert_allclose(beta, 1.0) + + +@pytest.mark.parametrize('axis', [0, 1, 2]) +@pytest.mark.parametrize('polarity', [-1, 1]) +def test_graded_expansion_uses_each_components_actual_yee_intervals(axis, polarity) -> None: + shape = [1, 1, 1] + shape[axis] = 4 + dxes = [[numpy.ones(n) for n in shape] for _ in range(2)] + dxes[0][axis] = numpy.array([1.0, 2.0, 3.0, 4.0]) + dxes[1][axis] = numpy.array([0.6, 0.8, 1.0, 1.2]) + source = [slice(None)] * 3 + source[axis] = slice(1, 2) + field = numpy.zeros((3, *shape), dtype=complex) + field[(slice(None), *source)] = 1 + 0.5j + expanded = waveguide_3d.expand_e(field, 0.2, dxes, axis, polarity, source) + for component in range(3): + line = expanded[component].ravel() + intervals = dxes[1][axis][1:] if component == axis else dxes[0][axis][:-1] + numpy.testing.assert_allclose(line[1:] / line[:-1], numpy.exp(-1j * polarity * 0.2 * intervals)) + numpy.testing.assert_array_equal(expanded[(slice(None), *source)], field[(slice(None), *source)]) + + +@pytest.mark.parametrize('beta', [2.0, 3.0]) +def test_grid_band_edge_rejects_power_normalization_and_selects_decay(beta) -> None: + shape = (1, 1, 5) + dxes = [[numpy.ones(n) for n in shape] for _ in range(2)] + kwargs = dict(omega=1.0, dxes=dxes, axis=2, polarity=1, + slices=(slice(None), slice(None), slice(0, 1)), epsilon=numpy.full((3, *shape), beta ** 2)) + with pytest.raises(ValueError, match='grid band edge'): + waveguide_3d.mode_from_2d(numpy.array([1, 0], dtype=complex), beta, **kwargs) + mode = waveguide_3d.mode_from_2d(numpy.array([1, 0], dtype=complex), beta, normalization='lorentz', **kwargs) + assert mode['wavenumber'].imag <= 0 + expanded = waveguide_3d.expand_e(mode['E'], mode['wavenumber'], dxes, 2, 1, kwargs['slices']) + ratio = abs(expanded[0, 0, 0, 1] / expanded[0, 0, 0, 0]) + assert ratio <= 1 + 1e-14 + if beta > 2: + assert ratio < 1 def test_waveguide_cyl_exy2e_and_exy2h_return_finite_full_fields() -> None: