7b01ddd789
Require a finite negative ln_R_per_layer before coefficient evaluation. Document None entries as the supported way to disable a face, and cover zero, signed zero, positive, and nonfinite targets with and without frequency shift. Validation: 209 focused tests passed; Ruff passed.
409 lines
20 KiB
Python
409 lines
20 KiB
Python
"""Independent regressions for the numerical/API feedback audit."""
|
|
import dataclasses
|
|
|
|
import numpy as np
|
|
import pytest
|
|
from numpy.testing import assert_allclose
|
|
from scipy import linalg, sparse
|
|
|
|
from .. import eigensolvers
|
|
from ..fdmath import vec
|
|
from ..fdmath.functional import deriv_forward, deriv_back
|
|
from ..fdtd import pml, base
|
|
from ..fdfd import eme, scpml, waveguide_2d as wg, waveguide_3d
|
|
|
|
|
|
def grid(shape):
|
|
return [[np.ones(n) for n in shape] for _ in range(2)]
|
|
|
|
|
|
def mode(profile, admittance=1, scale=1):
|
|
profile = np.asarray(profile)
|
|
e = np.zeros((3, profile.size, 1), complex)
|
|
h = np.zeros_like(e)
|
|
e[0, :, 0] = np.sqrt(2) * scale * profile
|
|
h[1, :, 0] = np.sqrt(2) * scale * admittance * profile
|
|
return vec(e), vec(h)
|
|
|
|
|
|
def test_disabled_cpml_matches_unequal_nonuniform_yee_updates():
|
|
rng = np.random.default_rng(6)
|
|
shape = (3, 4, 5, 6)
|
|
widths = [[np.linspace(.7 + .3 * family, 1.4 + .2 * family, n)
|
|
for n in shape[1:]] for family in range(2)]
|
|
e = rng.normal(size=shape) + 1j * rng.normal(size=shape)
|
|
h = rng.normal(size=shape) + 1j * rng.normal(size=shape)
|
|
eps, mu = 1 + rng.random(shape), 1 + rng.random(shape)
|
|
ec, hc = e.copy(), h.copy()
|
|
ue, uh = pml.updates_with_cpml([[None, None] for _ in range(3)], .1, widths, eps, dtype=complex)
|
|
ue(ec, hc, eps)
|
|
base.maxwell_e(.1, widths)(e, h, eps)
|
|
uh(ec, hc, mu)
|
|
base.maxwell_h(.1, widths)(e, h, mu)
|
|
assert_allclose(ec, e, atol=1e-14)
|
|
assert_allclose(hc, h, atol=1e-14)
|
|
|
|
|
|
def test_cpml_physical_spacing_scales_attenuation_rate():
|
|
small = pml.cpml_params(0, 1, dt=.1, dx=.25)
|
|
unit = pml.cpml_params(0, 1, dt=.4)
|
|
for family in ('param_e', 'param_h'):
|
|
assert_allclose(small[family], unit[family])
|
|
for invalid in (0, -1, np.inf, np.nan):
|
|
with pytest.raises(ValueError, match='dx'):
|
|
pml.cpml_params(0, 1, .1, dx=invalid)
|
|
|
|
|
|
def test_stretched_metric_satisfies_complex_reciprocity():
|
|
rng = np.random.default_rng(8)
|
|
shape = (4, 5)
|
|
widths = [[.7 + rng.random(n) - .2j * rng.random(n) for n in shape] for _ in range(2)]
|
|
operator = wg.operator_beta(.8, widths, vec(1 + rng.random((3, *shape))))
|
|
metric = wg.lorentz_metric(widths)
|
|
defect = operator.T @ metric - metric @ operator
|
|
assert sparse.linalg.norm(defect) / sparse.linalg.norm(metric @ operator) < 1e-12
|
|
physical = wg.lorentz_metric(widths, measure='physical')
|
|
assert sparse.linalg.norm(operator.T @ physical - physical @ operator) > 1e-3
|
|
e, h = mode([1, .3j, .2, .4])
|
|
# Direct one-dimensional complex-area quadrature, independent of the metric.
|
|
dx = [[np.full(4, 1 - .2j), np.ones(1)] for _ in range(2)]
|
|
assert_allclose(wg.inner_product(e, h, dx, measure='stretched'), (1 - .2j) * np.sum(np.array([1, .3j, .2, .4]) ** 2))
|
|
assert_allclose(wg.inner_product(e, h, dx, conj_h=True), np.sum(abs(np.array([1, .3j, .2, .4])) ** 2))
|
|
|
|
|
|
@pytest.mark.parametrize('magnetic', [False, True])
|
|
def test_reconstruction_can_normalize_complex_self_reaction(magnetic):
|
|
widths = [[np.full(2, 1 - .2j), np.ones(2)] for _ in range(2)]
|
|
transverse = np.concatenate((np.zeros(4), np.ones(4))) if magnetic else np.concatenate((np.ones(4), np.zeros(4)))
|
|
reconstruct = wg.normalized_fields_h if magnetic else wg.normalized_fields_e
|
|
e, h = reconstruct(transverse, 1, 1, widths, np.ones(12), normalization='lorentz')
|
|
assert_allclose(wg.inner_product(e, h, widths, measure='stretched'), 1, atol=1e-13)
|
|
with pytest.raises(ValueError, match='finite, nonzero'):
|
|
reconstruct(np.zeros(8), 1, 1, widths, np.ones(12))
|
|
with pytest.raises(ValueError, match='positive forward power'):
|
|
reconstruct(transverse, -1, 1, widths, np.ones(12))
|
|
|
|
|
|
@pytest.mark.parametrize('polarity', [-1, 1])
|
|
def test_scpml_face_copy_real_profile_and_final_stretch_cap(polarity):
|
|
widths = grid((10, 2))
|
|
result = scpml.stretch_with_scpml(
|
|
widths, 0, polarity, .2, thickness=3,
|
|
kappa_function=lambda depth: 1 + depth ** 2, max_imag_ratio=1,
|
|
)
|
|
selected = slice(0, 3) if polarity == -1 else slice(-3, None)
|
|
untouched = slice(3, None) if polarity == -1 else slice(None, -3)
|
|
assert result is not widths
|
|
for family in range(2):
|
|
assert_allclose(widths[family][0], 1)
|
|
assert_allclose(result[family][0][untouched], 1)
|
|
assert np.any(result[family][0][selected].real > 1)
|
|
assert np.any(result[family][0][selected].imag < 0)
|
|
assert np.all(abs(result[family][0].imag) <= result[family][0].real)
|
|
assert not np.shares_memory(result[family][1], widths[family][1])
|
|
mutated = scpml.stretch_with_scpml(widths, 0, polarity, .2, thickness=3, copy=False)
|
|
assert mutated is widths
|
|
assert np.any(widths[0][0].imag < 0)
|
|
|
|
|
|
def test_uniform_scpml_caps_final_dimensionless_stretch():
|
|
widths = scpml.uniform_grid_scpml((10, 1, 1), (3, 0, 0), .2,
|
|
kappa_function=lambda depth: 1 + depth, max_imag_ratio=.7)
|
|
for family in widths:
|
|
assert np.all(abs(family[0].imag) <= .7 * family[0].real + 1e-14)
|
|
assert_allclose(family[1], 1)
|
|
|
|
|
|
def test_indexed_eigensolve_accepts_seed_and_start_without_mutation():
|
|
operator = sparse.diags_array(np.linspace(-3, 5, 12))
|
|
first = eigensolvers.signed_eigensolve(operator, 2, rng=42)
|
|
second = eigensolvers.signed_eigensolve(operator, 2, rng=42)
|
|
assert_allclose(first[0], second[0], atol=1e-13)
|
|
assert_allclose(first[1], second[1], atol=1e-12)
|
|
start = np.linspace(1, 2, 12).astype(complex)
|
|
before = start.copy()
|
|
eigensolvers.signed_eigensolve(operator, 2, v0=start)
|
|
eigensolvers.power_iteration(operator, guess_vector=start)
|
|
assert_allclose(start, before)
|
|
with pytest.raises(ValueError, match='finite, nonzero'):
|
|
eigensolvers.signed_eigensolve(operator, 2, v0=np.zeros(12))
|
|
|
|
|
|
@pytest.mark.parametrize('axis', [0, 1, 2])
|
|
@pytest.mark.parametrize('polarity', [-1, 1])
|
|
def test_mode_from_2d_reuses_conversion_without_calling_eigensolver(monkeypatch, axis, polarity):
|
|
shape = (4, 5, 6)
|
|
widths, eps = grid(shape), np.ones((3, *shape))
|
|
slices = [slice(None)] * 3
|
|
slices[axis] = slice(1, 2)
|
|
points = int(np.prod([n for i, n in enumerate(shape) if i != axis]))
|
|
transverse = np.concatenate((np.ones(points), np.zeros(points)))
|
|
def forbidden(*args, **kwargs):
|
|
raise AssertionError('the conversion helper must not solve modes')
|
|
monkeypatch.setattr(wg, 'solve_mode', forbidden)
|
|
result = waveguide_3d.mode_from_2d(transverse, .5, omega=.5, dxes=widths, axis=axis,
|
|
polarity=polarity, slices=slices, epsilon=eps)
|
|
assert_allclose(result['wavenumber'], 2 * np.arcsin(.25))
|
|
assert np.isfinite(result['E']).all()
|
|
assert np.isfinite(result['H']).all()
|
|
monkeypatch.setattr(wg, 'solve_mode', lambda *args, **kwargs: (transverse, .5))
|
|
wrapped = waveguide_3d.solve_mode(0, .5, widths, axis, polarity, slices, eps)
|
|
assert_allclose(result['E'], wrapped['E'])
|
|
assert_allclose(result['H'], wrapped['H'])
|
|
|
|
|
|
def test_mixed_projection_scalar_overlap_and_unclipped_galerkin():
|
|
a = .8
|
|
left, right = mode([1, 0]), mode([a, np.sqrt(1 - a*a)])
|
|
ss = eme.get_s([left], [1], [right], [1], grid((2, 1)))
|
|
transmission, reflection = 2*a/(1+a*a), (1-a*a)/(1+a*a)
|
|
assert_allclose(ss, [[reflection, transmission], [transmission, -reflection]], atol=1e-14)
|
|
assert_allclose(ss.T @ ss, np.eye(2), atol=1e-14)
|
|
t, r = eme.get_tr([left], [1], [right], [1], grid((2, 1)), formulation='galerkin')
|
|
assert_allclose(t, 1/a)
|
|
assert_allclose(r, 0, atol=1e-14)
|
|
|
|
|
|
def test_asymmetric_transmission_uses_outgoing_rows_incoming_columns():
|
|
left, right = mode([1], scale=2), mode([1], scale=3)
|
|
ss = eme.get_s([left], [1], [right], [1], grid((1, 1)))
|
|
assert_allclose(ss, [[0, 1.5], [2/3, 0]], atol=1e-14)
|
|
# Continuity for the two physical incidence directions independently.
|
|
assert_allclose(2 * (1 + ss[0, 0]), 3 * ss[1, 0])
|
|
assert_allclose(3 * (1 + ss[1, 1]), 2 * ss[0, 1])
|
|
transfer = eme.get_abcd([left], [1], [right], [1], grid((1, 1))).toarray()
|
|
assert_allclose(transfer, np.eye(2) * (2/3), atol=1e-14)
|
|
assert_allclose(eme.abcd_to_s(transfer), ss, atol=1e-14)
|
|
assert_allclose(eme.abcd_to_s(sparse.csr_array(transfer)), ss, atol=1e-14)
|
|
|
|
|
|
def test_full_nondiagonal_gram_projections_match_exact_field_coordinates():
|
|
left_matrix = np.array([[1, .2j], [.1, 1]], complex)
|
|
right_matrix = np.array([[1, .3], [.2, 1]], complex)
|
|
dl = np.array([[1, .2], [.3j, 1]], complex)
|
|
dr = np.array([[1, -.1j], [.4, 1]], complex)
|
|
left = [mode(column) for column in left_matrix.T]
|
|
right = [mode(column) for column in right_matrix.T]
|
|
ss = eme.get_s(left, [1, 1], right, [1, 1], grid((2, 1)),
|
|
dual_ehLs=[mode(column) for column in dl.T],
|
|
dual_ehRs=[mode(column) for column in dr.T])
|
|
zero = np.zeros((2, 2))
|
|
expected = np.block([[zero, np.linalg.solve(left_matrix, right_matrix)],
|
|
[np.linalg.solve(right_matrix, left_matrix), zero]])
|
|
assert_allclose(ss, expected, atol=1e-13)
|
|
|
|
|
|
def test_mixed_interface_accepts_unequal_port_counts():
|
|
left = [mode([1, 0])]
|
|
right = [mode([1, 0]), mode([0, 1])]
|
|
ss = eme.get_s(left, [1], right, [1, 2], grid((2, 1)))
|
|
assert_allclose(ss, [[0, 1, 0], [1, 0, 0], [0, 0, -1]], atol=1e-14)
|
|
with pytest.raises(ValueError, match='equal port mode counts'):
|
|
eme.get_abcd(left, [1], right, [1, 2], grid((2, 1)))
|
|
|
|
|
|
def test_redheffer_matches_independent_internal_wave_elimination():
|
|
rng = np.random.default_rng(52)
|
|
left = .1 * (rng.normal(size=(3, 3)) + 1j * rng.normal(size=(3, 3)))
|
|
right = .1 * (rng.normal(size=(5, 5)) + 1j * rng.normal(size=(5, 5)))
|
|
internal = np.block([[np.eye(2), -left[1:, 1:]], [-right[:2, :2], np.eye(2)]])
|
|
incoming = np.block([[left[1:, :1], np.zeros((2, 3))], [np.zeros((2, 1)), right[:2, 2:]]])
|
|
waves = np.linalg.solve(internal, incoming)
|
|
outgoing = np.block([[left[:1, :1], np.zeros((1, 3))], [np.zeros((3, 1)), right[2:, 2:]]])
|
|
outgoing[:1] += left[:1, 1:] @ waves[2:]
|
|
outgoing[1:] += right[2:, :2] @ waves[:2]
|
|
assert_allclose(eme.cascade_s(left, right, internal_modes=2), outgoing, atol=1e-14)
|
|
|
|
|
|
def test_taper_overlap_has_no_generator_gain():
|
|
left, right = mode([1, 0]), mode([.8, .6])
|
|
sections = [eme.ModeSection(0, [left], [2]), eme.ModeSection(1, [right], [2])]
|
|
ss = eme.get_taper_s(sections, grid((2, 1)))
|
|
assert_allclose(np.linalg.svd(ss, compute_uv=False), 1, atol=1e-14)
|
|
assert_allclose(abs(ss[1, 0]), 2*.8/(1+.8**2), atol=1e-14)
|
|
with pytest.raises(TypeError, match='rtol'):
|
|
eme.get_taper_s(sections, grid((2, 1)), rtol=1e-8)
|
|
|
|
|
|
def test_long_leaky_taper_stays_finite_without_transfer_growth():
|
|
fields = mode([1])
|
|
sections = [eme.ModeSection(0, [fields], [1 - 1j]), eme.ModeSection(1000, [fields], [1 - 1j])]
|
|
ss = eme.get_taper_s(sections, grid((1, 1)))
|
|
assert_allclose(ss, 0, atol=1e-300)
|
|
with pytest.raises(np.linalg.LinAlgError, match='S to ABCD'):
|
|
eme.get_taper_abcd(sections, grid((1, 1)))
|
|
|
|
|
|
def test_prepared_dense_beta_overrides_are_aligned_and_avoid_overlaps(monkeypatch):
|
|
fields = [mode([1, 0]), mode([0, 1])]
|
|
phases = np.exp(1j * np.array([.3, -.7]))
|
|
shifted = [(e * phase, h * phase) for (e, h), phase in zip(fields, phases, strict=True)]
|
|
beta = np.array([[1 - .1j, .2], [.1j, 2 - .1j]])
|
|
sections = [eme.ModeSection(0, fields, [1, 2], wavenumber_matrix=beta),
|
|
eme.ModeSection(2, shifted, [1, 2], wavenumber_matrix=beta)]
|
|
prepared = eme.prepare_taper(sections, grid((2, 1)))
|
|
overrides = np.array([beta * 1.1, beta * .9])
|
|
changed = [dataclasses.replace(section, z=z, wavenumber_matrix=b)
|
|
for section, z, b in zip(sections, [0, 3], overrides, strict=True)]
|
|
expected = eme.get_taper_s(changed, grid((2, 1)))
|
|
def forbidden(*args, **kwargs):
|
|
raise AssertionError('evaluation must not recompute modal overlaps')
|
|
monkeypatch.setattr(eme, 'inner_product', forbidden)
|
|
monkeypatch.setattr(eme, 'lorentz_overlap', forbidden)
|
|
actual = eme.evaluate_taper_s(prepared, z_coords=[0, 3], wavenumber_matrices=overrides)
|
|
assert_allclose(actual, expected, atol=1e-13)
|
|
fields[0][0][:] = 100
|
|
assert_allclose(eme.evaluate_taper_s(prepared, z_coords=[0, 3], wavenumber_matrices=overrides), expected, atol=1e-13)
|
|
with pytest.raises(ValueError, match='read-only'):
|
|
prepared.interfaces[0, 0, 0] = 5
|
|
|
|
|
|
def test_explicit_backward_override_and_dense_uniform_propagation():
|
|
fields = [mode([1, 0]), mode([0, 1])]
|
|
backward = [(e, -h) for e, h in fields]
|
|
beta = np.array([[1 - .1j, .2], [.1j, 2 - .1j]])
|
|
minus = -beta + .03j * np.eye(2)
|
|
sections = [eme.ModeSection(z, fields, [1, 2], backward_modes=backward,
|
|
backward_wavenumbers=[-1, -2], wavenumber_matrix=beta,
|
|
backward_wavenumber_matrix=minus) for z in (0, 2)]
|
|
prepared = eme.prepare_taper(sections, grid((2, 1)))
|
|
with pytest.raises(ValueError, match='backward overrides'):
|
|
eme.evaluate_taper_s(prepared, wavenumber_matrices=np.array([beta, beta]))
|
|
ss = eme.evaluate_taper_s(prepared, wavenumber_matrices=np.array([beta, beta]),
|
|
backward_wavenumber_matrices=np.array([minus, minus]))
|
|
expected = np.block([[np.zeros((2, 2)), linalg.expm(2j * minus)],
|
|
[linalg.expm(-2j * beta), np.zeros((2, 2))]])
|
|
assert_allclose(ss, expected, atol=1e-13)
|
|
|
|
|
|
@pytest.mark.parametrize('dtype', [np.int8, np.float32, np.float64, np.complex64, np.complex128])
|
|
@pytest.mark.parametrize('axis', [0, 1, 2])
|
|
@pytest.mark.parametrize('forward', [False, True])
|
|
@pytest.mark.parametrize('spaced', [False, True])
|
|
@pytest.mark.parametrize('contiguous', [False, True])
|
|
def test_slice_derivatives_match_periodic_reference_and_output_buffers(dtype, axis, forward, spaced, contiguous):
|
|
f = np.arange(6 * 8 * 4).reshape(6, 8, 4).astype(dtype)[::2, ::2, ::2]
|
|
if np.issubdtype(dtype, np.complexfloating):
|
|
f = f * (1 + .3j)
|
|
if contiguous:
|
|
f = np.ascontiguousarray(f)
|
|
widths = [np.linspace(.7, 1.3, n) for n in f.shape] if spaced else None
|
|
derivative = (deriv_forward if forward else deriv_back)(widths)[axis]
|
|
expected = np.roll(f, -1, axis=axis) - f if forward else f - np.roll(f, 1, axis=axis)
|
|
if spaced:
|
|
shape = [1, 1, 1]
|
|
shape[axis] = f.shape[axis]
|
|
expected = expected / widths[axis].reshape(shape)
|
|
assert_allclose(derivative(f), expected, rtol=1e-6, atol=1e-6)
|
|
output = np.empty_like(expected)
|
|
assert derivative(f, out=output) is output
|
|
assert_allclose(output, expected, rtol=1e-6, atol=1e-6)
|
|
with pytest.raises(ValueError, match='alias'):
|
|
derivative(f, out=f)
|
|
singleton = np.ones((1, 1, 1), dtype=dtype)
|
|
plain = (deriv_forward if forward else deriv_back)()[axis]
|
|
assert_allclose(plain(singleton), 0)
|
|
|
|
|
|
def test_pml_cap_removes_above_light_line_artifact_in_homogeneous_fixture():
|
|
omega = .4
|
|
strong = grid((24, 1))
|
|
capped = grid((24, 1))
|
|
for polarity in (-1, 1):
|
|
strong = scpml.stretch_with_scpml(strong, 0, polarity, omega, thickness=5)
|
|
capped = scpml.stretch_with_scpml(capped, 0, polarity, omega, thickness=5, max_imag_ratio=1)
|
|
strong_values = np.linalg.eigvals(wg.operator_e(omega, strong, np.ones(72)).toarray())
|
|
capped_values = np.linalg.eigvals(wg.operator_e(omega, capped, np.ones(72)).toarray())
|
|
assert strong_values.real.max() > omega**2 + 1e-4
|
|
assert capped_values.real.max() <= omega**2 + 1e-10
|
|
|
|
|
|
def test_sampled_taper_converges_to_continuous_maxwell_propagation():
|
|
from scipy.integrate import solve_ivp
|
|
|
|
omega, length = 1.3, 2.0
|
|
def maxwell(z, state):
|
|
index = 1 + z / length
|
|
operator = np.array([[0, -1j * omega], [-1j * omega * index**2, 0]])
|
|
return (operator @ state.reshape(2, 2)).ravel()
|
|
integrated = solve_ivp(maxwell, (0, length), np.eye(2, dtype=complex).ravel(), rtol=1e-11, atol=1e-13)
|
|
left_fields = np.array([[1, 1], [1, -1]], complex)
|
|
right_fields = np.array([[1, 1], [2, -2]], complex) / np.sqrt(2)
|
|
transfer = np.linalg.solve(right_fields, integrated.y[:, -1].reshape(2, 2) @ left_fields)
|
|
a, b, c, d = transfer.ravel()
|
|
reference = np.array([[-c/d, 1/d], [a-b*c/d, b/d]])
|
|
errors = []
|
|
for count in (9, 17, 33):
|
|
sections = []
|
|
for z in np.linspace(0, length, count):
|
|
index = 1 + z / length
|
|
fields = mode([1], admittance=index, scale=1/np.sqrt(index))
|
|
sections.append(eme.ModeSection(z, [fields], [omega * index]))
|
|
scattering = eme.get_taper_s(sections, grid((1, 1)))
|
|
assert_allclose(scattering.conj().T @ scattering, np.eye(2), atol=1e-12)
|
|
errors.append(np.linalg.norm(scattering - reference))
|
|
assert errors[1] < errors[0] / 3
|
|
assert errors[2] < errors[1] / 3
|
|
assert errors[2] < 5e-4
|
|
|
|
|
|
def test_singular_conversions_and_interface_report_errors():
|
|
with pytest.raises(np.linalg.LinAlgError, match='S to ABCD'):
|
|
eme.s_to_abcd(np.eye(2))
|
|
with pytest.raises(np.linalg.LinAlgError, match='ABCD to S'):
|
|
eme.abcd_to_s(np.array([[1, 0], [0, 0]]))
|
|
fields = mode([0])
|
|
with pytest.raises(np.linalg.LinAlgError, match='interface continuity'):
|
|
eme.get_s([fields], [1], [fields], [1], grid((1, 1)))
|
|
mirror = np.eye(2)
|
|
with pytest.raises(np.linalg.LinAlgError, match='Redheffer internal feedback'):
|
|
eme.cascade_s(mirror, mirror, internal_modes=1)
|
|
|
|
|
|
def test_uniform_custom_profiles_leave_physical_aperture_unchanged():
|
|
widths = scpml.uniform_grid_scpml((10, 1, 1), (2, 0, 0), 1,
|
|
kappa_function=lambda depth: np.full_like(depth, 2))
|
|
for family in widths:
|
|
assert_allclose(family[0][3:8], 1)
|
|
|
|
|
|
@pytest.mark.parametrize('forward', [False, True])
|
|
def test_derivatives_do_not_subtract_across_independent_rows(forward):
|
|
fields = np.array([[[3e38, 3e38], [-3e38, -3e38]]], dtype=np.float32)
|
|
derivative = (deriv_forward if forward else deriv_back)()[2]
|
|
with np.errstate(over='raise', invalid='raise'):
|
|
assert_allclose(derivative(fields), 0)
|
|
|
|
|
|
@pytest.mark.parametrize(('shape', 'thicknesses'), [
|
|
((6, 8, 4), (0, 0, 0)),
|
|
((6, 8, 4), (1, 1, 1)),
|
|
((6, 8, 4), (2, 3, 0)),
|
|
((6, 8, 4), (3, 4, 2)),
|
|
])
|
|
@pytest.mark.parametrize('kappa', [None, 'polynomial', 'constant'])
|
|
@pytest.mark.parametrize('cap', [None, .7, 1.0])
|
|
def test_scpml_builders_agree_including_real_profiles_and_touching_faces(shape, thicknesses, kappa, cap):
|
|
profile = scpml.prepare_s_function(ln_R=-9, m=2.5)
|
|
real_profile = {
|
|
None: None,
|
|
'polynomial': lambda depth: 1 + 2 * depth**2,
|
|
'constant': lambda depth: np.full_like(depth, 2),
|
|
}[kappa]
|
|
options = dict(omega=.8, epsilon_effective=2.25, s_function=profile,
|
|
kappa_function=real_profile, max_imag_ratio=cap)
|
|
uniform = scpml.uniform_grid_scpml(shape, thicknesses, **options)
|
|
explicit = grid(shape)
|
|
for axis, thickness in enumerate(thicknesses):
|
|
for polarity in (-1, 1):
|
|
explicit = scpml.stretch_with_scpml(explicit, axis, polarity, thickness=thickness, **options)
|
|
for family in range(2):
|
|
for axis in range(3):
|
|
assert_allclose(uniform[family][axis], explicit[family][axis], rtol=1e-14, atol=1e-14)
|
|
|
|
|
|
@pytest.mark.parametrize('attenuation', [0.0, -0.0, .1, np.inf, -np.inf, np.nan])
|
|
@pytest.mark.parametrize('alpha', [0, .2])
|
|
def test_cpml_rejects_zero_or_invalid_attenuation_before_coefficient_evaluation(attenuation, alpha):
|
|
with np.errstate(divide='raise', invalid='raise'), pytest.raises(ValueError, match='ln_R_per_layer must be finite and negative'):
|
|
pml.cpml_params(0, 1, .1, ln_R_per_layer=attenuation, cfs_alpha=alpha)
|