jan/meanas
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

[tests] refactor tests

Browse source
This commit is contained in:
Jan Petykiewicz 2026-04-18 00:52:04 -07:00
commit 8cdcd08ba0
25 changed files with 649 additions and 616 deletions
Download patch file Download diff file Expand all files Collapse all files

37
meanas/test/_bloch_case.py Normal file
View file

@ -0,0 +1,37 @@
import numpy
SHAPE = (2, 2, 2)
G_MATRIX = numpy.eye(3)
K0_GENERAL = numpy.array([0.1, 0.2, 0.3], dtype=float)
K0_AXIAL = numpy.array([0.0, 0.0, 0.25], dtype=float)
K0_X = numpy.array([0.1, 0.0, 0.0], dtype=float)
EPSILON = numpy.ones((3, *SHAPE), dtype=float)
MU = numpy.stack([
numpy.linspace(2.0, 2.7, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(2.1, 2.8, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(2.2, 2.9, numpy.prod(SHAPE)).reshape(SHAPE),
])
H_SIZE = 2 * numpy.prod(SHAPE)
H_MN = (numpy.arange(H_SIZE) + 0.25j).astype(complex)
ZERO_H_MN = numpy.zeros_like(H_MN)
Y0 = (numpy.arange(H_SIZE, dtype=float) + 1j * numpy.linspace(0.1, 0.9, H_SIZE))[None, :]
Y0_TWO_MODE = numpy.vstack(
[
numpy.arange(H_SIZE, dtype=float) + 1j * numpy.linspace(0.1, 0.9, H_SIZE),
numpy.linspace(2.0, 3.5, H_SIZE) - 0.5j * numpy.arange(H_SIZE, dtype=float),
],
)
def build_overlap_fixture() -> tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]:
e_in = numpy.zeros((3, *SHAPE), dtype=complex)
h_in = numpy.zeros_like(e_in)
e_out = numpy.zeros_like(e_in)
h_out = numpy.zeros_like(e_in)
e_in[1] = 1.0
h_in[2] = 2.0
e_out[1] = 3.0
h_out[2] = 4.0
return e_in, h_in, e_out, h_out

49
meanas/test/_fdfd_case.py Normal file
View file

@ -0,0 +1,49 @@
import numpy
from ..fdmath import vec, unvec
OMEGA = 1 / 1500
SHAPE = (2, 3, 2)
DXES = [
[numpy.array([1.0, 1.5]), numpy.array([0.75, 1.25, 1.5]), numpy.array([1.2, 0.8])],
[numpy.array([0.9, 1.4]), numpy.array([0.8, 1.1, 1.4]), numpy.array([1.0, 0.7])],
]
EPSILON = numpy.stack([
numpy.linspace(1.0, 2.2, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(1.1, 2.3, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(1.2, 2.4, numpy.prod(SHAPE)).reshape(SHAPE),
])
MU = numpy.stack([
numpy.linspace(2.0, 3.2, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(2.1, 3.3, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(2.2, 3.4, numpy.prod(SHAPE)).reshape(SHAPE),
])
E_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) + 0.5j).astype(complex)
H_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) * 0.25 - 0.75j).astype(complex)
PEC = numpy.zeros((3, *SHAPE), dtype=float)
PEC[1, 0, 1, 0] = 1.0
PMC = numpy.zeros((3, *SHAPE), dtype=float)
PMC[2, 1, 2, 1] = 1.0
TF_REGION = numpy.zeros((3, *SHAPE), dtype=float)
TF_REGION[:, 0, 1, 0] = 1.0
BOUNDARY_SHAPE = (3, 4, 3)
BOUNDARY_DXES = [
[numpy.array([1.0, 1.5, 0.8]), numpy.array([0.75, 1.25, 1.5, 0.9]), numpy.array([1.2, 0.8, 1.1])],
[numpy.array([0.9, 1.4, 1.0]), numpy.array([0.8, 1.1, 1.4, 1.0]), numpy.array([1.0, 0.7, 1.3])],
]
BOUNDARY_EPSILON = numpy.stack([
numpy.linspace(1.0, 2.2, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
numpy.linspace(1.1, 2.3, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
numpy.linspace(1.2, 2.4, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
])
BOUNDARY_FIELD = (numpy.arange(3 * numpy.prod(BOUNDARY_SHAPE)).reshape((3, *BOUNDARY_SHAPE)) + 0.5j).astype(complex)
def apply_fdfd_matrix(op, field: numpy.ndarray, shape: tuple[int, ...] = SHAPE) -> numpy.ndarray:
return unvec(op @ vec(field), shape)

56
meanas/test/_solver_cases.py Normal file
View file

@ -0,0 +1,56 @@
import dataclasses
import numpy
from scipy import sparse
import scipy.sparse.linalg as spalg
from ._test_builders import unit_dxes
@dataclasses.dataclass(frozen=True)
class DiagonalEigenCase:
operator: sparse.csr_matrix
linear_operator: spalg.LinearOperator
guess_default: numpy.ndarray
guess_sparse: numpy.ndarray
@dataclasses.dataclass(frozen=True)
class SolverPlumbingCase:
omega: float
a0: sparse.csr_matrix
pl: sparse.csr_matrix
pr: sparse.csr_matrix
j: numpy.ndarray
guess: numpy.ndarray
solver_result: numpy.ndarray
dxes: tuple[tuple[numpy.ndarray, ...], tuple[numpy.ndarray, ...]]
epsilon: numpy.ndarray
def diagonal_eigen_case() -> DiagonalEigenCase:
operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
linear_operator = spalg.LinearOperator(
shape=operator.shape,
dtype=complex,
matvec=lambda vv: operator @ vv,
)
return DiagonalEigenCase(
operator=operator,
linear_operator=linear_operator,
guess_default=numpy.array([0.0, 1.0, 1e-6, 0.0], dtype=complex),
guess_sparse=numpy.array([1.0, 0.1, 0.0, 0.0], dtype=complex),
)
def solver_plumbing_case() -> SolverPlumbingCase:
return SolverPlumbingCase(
omega=2.0,
a0=sparse.csr_matrix(numpy.array([[1.0 + 2.0j, 2.0], [3.0 - 1.0j, 4.0]])),
pl=sparse.csr_matrix(numpy.array([[2.0, 0.0], [0.0, 3.0j]])),
pr=sparse.csr_matrix(numpy.array([[0.5, 0.0], [0.0, -2.0j]])),
j=numpy.array([1.0 + 0.5j, -2.0]),
guess=numpy.array([0.25 - 0.75j, 1.5 + 0.5j]),
solver_result=numpy.array([3.0 - 1.0j, -4.0 + 2.0j]),
dxes=unit_dxes((1, 1, 1)),
epsilon=numpy.ones(2),
)

22
meanas/test/_test_builders.py Normal file
View file

@ -0,0 +1,22 @@
import numpy
def real_ramp(shape: tuple[int, ...], *, scale: float = 1.0, offset: float = 0.0) -> numpy.ndarray:
return numpy.arange(numpy.prod(shape), dtype=float).reshape(shape, order='C') * scale + offset
def complex_ramp(
shape: tuple[int, ...],
*,
scale: float = 1.0,
offset: float = 0.0,
imag_scale: float = 0.0,
imag_offset: float = 0.0,
) -> numpy.ndarray:
real = real_ramp(shape, scale=scale, offset=offset)
imag = real_ramp(shape, scale=imag_scale, offset=imag_offset)
return (real + 1j * imag).astype(complex)
def unit_dxes(shape: tuple[int, ...]) -> tuple[tuple[numpy.ndarray, ...], tuple[numpy.ndarray, ...]]:
return tuple(tuple(numpy.ones(length) for length in shape) for _ in range(2)) # type: ignore[return-value]

15
meanas/test/conftest.py
View file

@ -9,7 +9,7 @@ import numpy
from numpy.typing import NDArray from numpy.typing import NDArray
import pytest # type: ignore import pytest # type: ignore
from .utils import PRNG from .utils import make_prng
FixtureRequest = Any FixtureRequest = Any
@ -42,6 +42,7 @@ def epsilon(
epsilon_bg: float, epsilon_bg: float,
epsilon_fg: float, epsilon_fg: float,
) -> NDArray[numpy.float64]: ) -> NDArray[numpy.float64]:
prng = make_prng()
is3d = (numpy.array(shape) == 1).sum() == 0 is3d = (numpy.array(shape) == 1).sum() == 0
if is3d: if is3d:
if request.param == '000': if request.param == '000':
@ -57,9 +58,11 @@ def epsilon(
elif request.param == '000': elif request.param == '000':
epsilon[:, 0, 0, 0] = epsilon_fg epsilon[:, 0, 0, 0] = epsilon_fg
elif request.param == 'random': elif request.param == 'random':
epsilon[:] = PRNG.uniform(low=min(epsilon_bg, epsilon_fg), epsilon[:] = prng.uniform(
high=max(epsilon_bg, epsilon_fg), low=min(epsilon_bg, epsilon_fg),
size=shape) high=max(epsilon_bg, epsilon_fg),
size=shape,
)
return epsilon return epsilon
@ -80,6 +83,7 @@ def dxes(
shape: tuple[int, ...], shape: tuple[int, ...],
dx: float, dx: float,
) -> list[list[NDArray[numpy.float64]]]: ) -> list[list[NDArray[numpy.float64]]]:
prng = make_prng()
if request.param == 'uniform': if request.param == 'uniform':
dxes = [[numpy.full(s, dx) for s in shape[1:]] for _ in range(2)] dxes = [[numpy.full(s, dx) for s in shape[1:]] for _ in range(2)]
elif request.param == 'centerbig': elif request.param == 'centerbig':
@ -88,8 +92,7 @@ def dxes(
for ax in (0, 1, 2): for ax in (0, 1, 2):
dxes[eh][ax][dxes[eh][ax].size // 2] *= 1.1 dxes[eh][ax][dxes[eh][ax].size // 2] *= 1.1
elif request.param == 'random': elif request.param == 'random':
dxe = [PRNG.uniform(low=1.0 * dx, high=1.1 * dx, size=s) for s in shape[1:]] dxe = [prng.uniform(low=1.0 * dx, high=1.1 * dx, size=s) for s in shape[1:]]
dxh = [(d + numpy.roll(d, -1)) / 2 for d in dxe] dxh = [(d + numpy.roll(d, -1)) / 2 for d in dxe]
dxes = [dxe, dxh] dxes = [dxe, dxh]
return dxes return dxes

48
meanas/test/test_bloch_foundations.py
View file

@ -2,23 +2,11 @@ import numpy
from numpy.linalg import norm from numpy.linalg import norm
from ..fdfd import bloch from ..fdfd import bloch
from ._bloch_case import EPSILON, G_MATRIX, H_MN, K0_AXIAL, K0_GENERAL, MU, SHAPE, ZERO_H_MN
from .utils import assert_close
SHAPE = (2, 2, 2)
K0 = numpy.array([0.1, 0.2, 0.3])
G_MATRIX = numpy.eye(3)
EPSILON = numpy.ones((3, *SHAPE), dtype=float)
MU = numpy.stack([
numpy.linspace(2.0, 2.7, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(2.1, 2.8, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(2.2, 2.9, numpy.prod(SHAPE)).reshape(SHAPE),
])
H_MN = (numpy.arange(2 * numpy.prod(SHAPE)) + 0.25j).astype(complex)
ZERO_H_MN = numpy.zeros_like(H_MN)
def test_generate_kmn_general_case_returns_orthonormal_basis() -> None: def test_generate_kmn_general_case_returns_orthonormal_basis() -> None:
k_mag, m_vecs, n_vecs = bloch.generate_kmn(K0, G_MATRIX, SHAPE) k_mag, m_vecs, n_vecs = bloch.generate_kmn(K0_GENERAL, G_MATRIX, SHAPE)
assert k_mag.shape == SHAPE + (1,) assert k_mag.shape == SHAPE + (1,)
assert m_vecs.shape == SHAPE + (3,) assert m_vecs.shape == SHAPE + (3,)
@ -27,57 +15,57 @@ def test_generate_kmn_general_case_returns_orthonormal_basis() -> None:
assert numpy.isfinite(m_vecs).all() assert numpy.isfinite(m_vecs).all()
assert numpy.isfinite(n_vecs).all() assert numpy.isfinite(n_vecs).all()
numpy.testing.assert_allclose(norm(m_vecs.reshape(-1, 3), axis=1), 1.0) assert_close(norm(m_vecs.reshape(-1, 3), axis=1), 1.0)
numpy.testing.assert_allclose(norm(n_vecs.reshape(-1, 3), axis=1), 1.0) assert_close(norm(n_vecs.reshape(-1, 3), axis=1), 1.0)
numpy.testing.assert_allclose(numpy.sum(m_vecs * n_vecs, axis=3), 0.0, atol=1e-12) assert_close(numpy.sum(m_vecs * n_vecs, axis=3), 0.0, atol=1e-12)
def test_generate_kmn_z_aligned_uses_default_transverse_basis() -> None: def test_generate_kmn_z_aligned_uses_default_transverse_basis() -> None:
k_mag, m_vecs, n_vecs = bloch.generate_kmn([0.0, 0.0, 0.25], G_MATRIX, (1, 1, 1)) k_mag, m_vecs, n_vecs = bloch.generate_kmn(K0_AXIAL, G_MATRIX, (1, 1, 1))
assert numpy.isfinite(k_mag).all() assert numpy.isfinite(k_mag).all()
numpy.testing.assert_allclose(m_vecs[0, 0, 0], [0.0, 1.0, 0.0]) assert_close(m_vecs[0, 0, 0], [0.0, 1.0, 0.0])
numpy.testing.assert_allclose(numpy.sum(m_vecs * n_vecs, axis=3), 0.0, atol=1e-12) assert_close(numpy.sum(m_vecs * n_vecs, axis=3), 0.0, atol=1e-12)
numpy.testing.assert_allclose(norm(n_vecs.reshape(-1, 3), axis=1), 1.0) assert_close(norm(n_vecs.reshape(-1, 3), axis=1), 1.0)
def test_maxwell_operator_returns_finite_column_vector_without_mu() -> None: def test_maxwell_operator_returns_finite_column_vector_without_mu() -> None:
operator = bloch.maxwell_operator(K0, G_MATRIX, EPSILON) operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON)
result = operator(H_MN.copy()) result = operator(H_MN.copy())
zero_result = operator(ZERO_H_MN.copy()) zero_result = operator(ZERO_H_MN.copy())
assert result.shape == (2 * numpy.prod(SHAPE), 1) assert result.shape == (2 * numpy.prod(SHAPE), 1)
assert numpy.isfinite(result).all() assert numpy.isfinite(result).all()
numpy.testing.assert_allclose(zero_result, 0.0) assert_close(zero_result, 0.0)
def test_maxwell_operator_returns_finite_column_vector_with_mu() -> None: def test_maxwell_operator_returns_finite_column_vector_with_mu() -> None:
operator = bloch.maxwell_operator(K0, G_MATRIX, EPSILON, MU) operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON, MU)
result = operator(H_MN.copy()) result = operator(H_MN.copy())
zero_result = operator(ZERO_H_MN.copy()) zero_result = operator(ZERO_H_MN.copy())
assert result.shape == (2 * numpy.prod(SHAPE), 1) assert result.shape == (2 * numpy.prod(SHAPE), 1)
assert numpy.isfinite(result).all() assert numpy.isfinite(result).all()
numpy.testing.assert_allclose(zero_result, 0.0) assert_close(zero_result, 0.0)
def test_inverse_maxwell_operator_returns_finite_column_vector_for_both_mu_branches() -> None: def test_inverse_maxwell_operator_returns_finite_column_vector_for_both_mu_branches() -> None:
for mu in (None, MU): for mu in (None, MU):
operator = bloch.inverse_maxwell_operator_approx(K0, G_MATRIX, EPSILON, mu) operator = bloch.inverse_maxwell_operator_approx(K0_GENERAL, G_MATRIX, EPSILON, mu)
result = operator(H_MN.copy()) result = operator(H_MN.copy())
zero_result = operator(ZERO_H_MN.copy()) zero_result = operator(ZERO_H_MN.copy())
assert result.shape == (2 * numpy.prod(SHAPE), 1) assert result.shape == (2 * numpy.prod(SHAPE), 1)
assert numpy.isfinite(result).all() assert numpy.isfinite(result).all()
numpy.testing.assert_allclose(zero_result, 0.0) assert_close(zero_result, 0.0)
def test_bloch_field_converters_return_finite_fields() -> None: def test_bloch_field_converters_return_finite_fields() -> None:
e_field = bloch.hmn_2_exyz(K0, G_MATRIX, EPSILON)(H_MN.copy()) e_field = bloch.hmn_2_exyz(K0_GENERAL, G_MATRIX, EPSILON)(H_MN.copy())
h_field = bloch.hmn_2_hxyz(K0, G_MATRIX, EPSILON)(H_MN.copy()) h_field = bloch.hmn_2_hxyz(K0_GENERAL, G_MATRIX, EPSILON)(H_MN.copy())
assert e_field.shape == (3, *SHAPE) assert e_field.shape == (3, *SHAPE)
assert h_field.shape == (3, *SHAPE) assert h_field.shape == (3, *SHAPE)

67
meanas/test/test_bloch_interactions.py
View file

@ -4,33 +4,8 @@ from numpy.testing import assert_allclose
from types import SimpleNamespace from types import SimpleNamespace
from ..fdfd import bloch from ..fdfd import bloch
from ._bloch_case import EPSILON, G_MATRIX, H_SIZE, K0_X, SHAPE, Y0, Y0_TWO_MODE, build_overlap_fixture
from .utils import assert_close
SHAPE = (2, 2, 2)
G_MATRIX = numpy.eye(3)
EPSILON = numpy.ones((3, *SHAPE), dtype=float)
K0 = numpy.array([0.1, 0.0, 0.0], dtype=float)
H_SIZE = 2 * numpy.prod(SHAPE)
Y0 = (numpy.arange(H_SIZE, dtype=float) + 1j * numpy.linspace(0.1, 0.9, H_SIZE))[None, :]
Y0_TWO_MODE = numpy.vstack(
[
numpy.arange(H_SIZE, dtype=float) + 1j * numpy.linspace(0.1, 0.9, H_SIZE),
numpy.linspace(2.0, 3.5, H_SIZE) - 0.5j * numpy.arange(H_SIZE, dtype=float),
],
)
def build_overlap_fixture() -> tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]:
e_in = numpy.zeros((3, *SHAPE), dtype=complex)
h_in = numpy.zeros_like(e_in)
e_out = numpy.zeros_like(e_in)
h_out = numpy.zeros_like(e_in)
e_in[1] = 1.0
h_in[2] = 2.0
e_out[1] = 3.0
h_out[2] = 4.0
return e_in, h_in, e_out, h_out
def test_rtrace_atb_matches_real_frobenius_inner_product() -> None: def test_rtrace_atb_matches_real_frobenius_inner_product() -> None:
@ -45,8 +20,8 @@ def test_symmetrize_returns_hermitian_average() -> None:
matrix = numpy.array([[1.0 + 2.0j, 3.0 - 1.0j], [2.0j, 4.0]], dtype=complex) matrix = numpy.array([[1.0 + 2.0j, 3.0 - 1.0j], [2.0j, 4.0]], dtype=complex)
result = bloch._symmetrize(matrix) result = bloch._symmetrize(matrix)
assert_allclose(result, 0.5 * (matrix + matrix.conj().T)) assert_close(result, 0.5 * (matrix + matrix.conj().T))
assert_allclose(result, result.conj().T) assert_close(result, result.conj().T)
def test_inner_product_is_nonmutating_and_obeys_sign_symmetry() -> None: def test_inner_product_is_nonmutating_and_obeys_sign_symmetry() -> None:
@ -58,9 +33,9 @@ def test_inner_product_is_nonmutating_and_obeys_sign_symmetry() -> None:
np_term = bloch.inner_product(e_out, -h_out, e_in, h_in) np_term = bloch.inner_product(e_out, -h_out, e_in, h_in)
nn = bloch.inner_product(e_out, -h_out, e_in, -h_in) nn = bloch.inner_product(e_out, -h_out, e_in, -h_in)
assert_allclose(pp, 0.8164965809277263 + 0.0j) assert_close(pp, 0.8164965809277263 + 0.0j)
assert_allclose(pp, -nn, atol=1e-12, rtol=1e-12) assert_close(pp, -nn, atol=1e-12, rtol=1e-12)
assert_allclose(pn, -np_term, atol=1e-12, rtol=1e-12) assert_close(pn, -np_term, atol=1e-12, rtol=1e-12)
assert numpy.array_equal(e_in, originals[0]) assert numpy.array_equal(e_in, originals[0])
assert numpy.array_equal(h_in, originals[1]) assert numpy.array_equal(h_in, originals[1])
assert numpy.array_equal(e_out, originals[2]) assert numpy.array_equal(e_out, originals[2])
@ -72,8 +47,8 @@ def test_trq_returns_expected_transmission_and_reflection() -> None:
transmission, reflection = bloch.trq(e_in, h_in, e_out, h_out) transmission, reflection = bloch.trq(e_in, h_in, e_out, h_out)
assert_allclose(transmission, 0.9797958971132713 + 0.0j, atol=1e-12, rtol=1e-12) assert_close(transmission, 0.9797958971132713 + 0.0j, atol=1e-12, rtol=1e-12)
assert_allclose(reflection, 0.2 + 0.0j, atol=1e-12, rtol=1e-12) assert_close(reflection, 0.2 + 0.0j, atol=1e-12, rtol=1e-12)
def test_eigsolve_returns_finite_modes_with_small_residual() -> None: def test_eigsolve_returns_finite_modes_with_small_residual() -> None:
@ -85,7 +60,7 @@ def test_eigsolve_returns_finite_modes_with_small_residual() -> None:
eigvals, eigvecs = bloch.eigsolve( eigvals, eigvecs = bloch.eigsolve(
1, 1,
K0, K0_X,
G_MATRIX, G_MATRIX,
EPSILON, EPSILON,
tolerance=1e-6, tolerance=1e-6,
@ -94,7 +69,7 @@ def test_eigsolve_returns_finite_modes_with_small_residual() -> None:
callback=callback, callback=callback,
) )
operator = bloch.maxwell_operator(K0, G_MATRIX, EPSILON) operator = bloch.maxwell_operator(K0_X, G_MATRIX, EPSILON)
eigvec = eigvecs[0] / numpy.linalg.norm(eigvecs[0]) eigvec = eigvecs[0] / numpy.linalg.norm(eigvecs[0])
residual = numpy.linalg.norm(operator(eigvec).reshape(-1) - eigvals[0] * eigvec) / numpy.linalg.norm(eigvec) residual = numpy.linalg.norm(operator(eigvec).reshape(-1) - eigvals[0] * eigvec) / numpy.linalg.norm(eigvec)
@ -109,7 +84,7 @@ def test_eigsolve_returns_finite_modes_with_small_residual() -> None:
def test_eigsolve_without_initial_guess_returns_finite_modes() -> None: def test_eigsolve_without_initial_guess_returns_finite_modes() -> None:
eigvals, eigvecs = bloch.eigsolve( eigvals, eigvecs = bloch.eigsolve(
1, 1,
K0, K0_X,
G_MATRIX, G_MATRIX,
EPSILON, EPSILON,
tolerance=1e-6, tolerance=1e-6,
@ -117,7 +92,7 @@ def test_eigsolve_without_initial_guess_returns_finite_modes() -> None:
y0=None, y0=None,
) )
operator = bloch.maxwell_operator(K0, G_MATRIX, EPSILON) operator = bloch.maxwell_operator(K0_X, G_MATRIX, EPSILON)
eigvec = eigvecs[0] / numpy.linalg.norm(eigvecs[0]) eigvec = eigvecs[0] / numpy.linalg.norm(eigvecs[0])
residual = numpy.linalg.norm(operator(eigvec).reshape(-1) - eigvals[0] * eigvec) / numpy.linalg.norm(eigvec) residual = numpy.linalg.norm(operator(eigvec).reshape(-1) - eigvals[0] * eigvec) / numpy.linalg.norm(eigvec)
@ -143,7 +118,7 @@ def test_eigsolve_recovers_from_singular_initial_subspace(monkeypatch: pytest.Mo
eigvals, eigvecs = bloch.eigsolve( eigvals, eigvecs = bloch.eigsolve(
2, 2,
K0, K0_X,
G_MATRIX, G_MATRIX,
EPSILON, EPSILON,
tolerance=1e-6, tolerance=1e-6,
@ -181,7 +156,7 @@ def test_eigsolve_reconditions_large_trace_initial_subspace(monkeypatch: pytest.
eigvals, eigvecs = bloch.eigsolve( eigvals, eigvecs = bloch.eigsolve(
2, 2,
K0, K0_X,
G_MATRIX, G_MATRIX,
EPSILON, EPSILON,
tolerance=1e-6, tolerance=1e-6,
@ -208,7 +183,7 @@ def test_eigsolve_qi_memoization_reuses_cached_theta(monkeypatch: pytest.MonkeyP
eigvals, eigvecs = bloch.eigsolve( eigvals, eigvecs = bloch.eigsolve(
1, 1,
K0, K0_X,
G_MATRIX, G_MATRIX,
EPSILON, EPSILON,
tolerance=1e-6, tolerance=1e-6,
@ -243,7 +218,7 @@ def test_eigsolve_qi_taylor_expansions_return_finite_modes(monkeypatch: pytest.M
eigvals, eigvecs = bloch.eigsolve( eigvals, eigvecs = bloch.eigsolve(
1, 1,
K0, K0_X,
G_MATRIX, G_MATRIX,
EPSILON, EPSILON,
tolerance=1e-6, tolerance=1e-6,
@ -278,7 +253,7 @@ def test_eigsolve_qi_inexplicable_singularity_raises(monkeypatch: pytest.MonkeyP
with pytest.raises(Exception, match='Inexplicable singularity in trace_func'): with pytest.raises(Exception, match='Inexplicable singularity in trace_func'):
bloch.eigsolve( bloch.eigsolve(
1, 1,
K0, K0_X,
G_MATRIX, G_MATRIX,
EPSILON, EPSILON,
tolerance=1e-6, tolerance=1e-6,
@ -290,7 +265,7 @@ def test_eigsolve_qi_inexplicable_singularity_raises(monkeypatch: pytest.MonkeyP
def test_find_k_returns_vector_frequency_and_callbacks() -> None: def test_find_k_returns_vector_frequency_and_callbacks() -> None:
target_eigvals, _target_eigvecs = bloch.eigsolve( target_eigvals, _target_eigvecs = bloch.eigsolve(
1, 1,
K0, K0_X,
G_MATRIX, G_MATRIX,
EPSILON, EPSILON,
tolerance=1e-6, tolerance=1e-6,
@ -329,8 +304,8 @@ def test_find_k_returns_vector_frequency_and_callbacks() -> None:
assert found_k.shape == (3,) assert found_k.shape == (3,)
assert numpy.isfinite(found_k).all() assert numpy.isfinite(found_k).all()
assert_allclose(numpy.cross(found_k, [1.0, 0.0, 0.0]), 0.0, atol=1e-12, rtol=1e-12) assert_close(numpy.cross(found_k, [1.0, 0.0, 0.0]), 0.0, atol=1e-12, rtol=1e-12)
assert_allclose(found_k, K0, atol=1e-4, rtol=1e-4) assert_close(found_k, K0_X, atol=1e-4, rtol=1e-4)
assert abs(found_frequency - target_frequency) <= 1e-4 assert abs(found_frequency - target_frequency) <= 1e-4
assert eigvals.shape == (1,) assert eigvals.shape == (1,)
assert eigvecs.shape == (1, H_SIZE) assert eigvecs.shape == (1, H_SIZE)

90
meanas/test/test_eigensolvers_numerics.py
View file

@ -1,46 +1,40 @@
import numpy import numpy
from numpy.linalg import norm from numpy.linalg import norm
from scipy import sparse import pytest
import scipy.sparse.linalg as spalg
from numpy.testing import assert_allclose
from ..eigensolvers import rayleigh_quotient_iteration, signed_eigensolve from ._solver_cases import diagonal_eigen_case
from .utils import assert_close
from ..eigensolvers import power_iteration, rayleigh_quotient_iteration, signed_eigensolve
def test_rayleigh_quotient_iteration_with_linear_operator() -> None: def test_rayleigh_quotient_iteration_with_linear_operator() -> None:
operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr() case = diagonal_eigen_case()
linear_operator = spalg.LinearOperator(
shape=operator.shape,
dtype=complex,
matvec=lambda vv: operator @ vv,
)
def dense_solver( def dense_solver(
shifted_operator: spalg.LinearOperator, shifted_operator,
rhs: numpy.ndarray, rhs: numpy.ndarray,
) -> numpy.ndarray: ) -> numpy.ndarray:
basis = numpy.eye(operator.shape[0], dtype=complex) basis = numpy.eye(case.operator.shape[0], dtype=complex)
columns = [shifted_operator.matvec(basis[:, ii]) for ii in range(operator.shape[0])] columns = [shifted_operator.matvec(basis[:, ii]) for ii in range(case.operator.shape[0])]
dense_matrix = numpy.column_stack(columns) dense_matrix = numpy.column_stack(columns)
return numpy.linalg.lstsq(dense_matrix, rhs, rcond=None)[0] return numpy.linalg.lstsq(dense_matrix, rhs, rcond=None)[0]
guess = numpy.array([0.0, 1.0, 1e-6, 0.0], dtype=complex)
eigval, eigvec = rayleigh_quotient_iteration( eigval, eigvec = rayleigh_quotient_iteration(
linear_operator, case.linear_operator,
guess, case.guess_default,
iterations=8, iterations=8,
solver=dense_solver, solver=dense_solver,
) )
residual = norm(operator @ eigvec - eigval * eigvec) residual = norm(case.operator @ eigvec - eigval * eigvec)
assert abs(eigval - 3.0) < 1e-12 assert abs(eigval - 3.0) < 1e-12
assert residual < 1e-12 assert residual < 1e-12
def test_signed_eigensolve_negative_returns_largest_negative_mode() -> None: def test_signed_eigensolve_negative_returns_largest_negative_mode() -> None:
operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr() case = diagonal_eigen_case()
eigvals, eigvecs = signed_eigensolve(operator, how_many=1, negative=True) eigvals, eigvecs = signed_eigensolve(case.operator, how_many=1, negative=True)
assert eigvals.shape == (1,) assert eigvals.shape == (1,)
assert eigvecs.shape == (4, 1) assert eigvecs.shape == (4, 1)
@ -49,35 +43,59 @@ def test_signed_eigensolve_negative_returns_largest_negative_mode() -> None:
def test_rayleigh_quotient_iteration_uses_default_linear_operator_solver() -> None: def test_rayleigh_quotient_iteration_uses_default_linear_operator_solver() -> None:
operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr() case = diagonal_eigen_case()
linear_operator = spalg.LinearOperator(
shape=operator.shape,
dtype=complex,
matvec=lambda vv: operator @ vv,
)
eigval, eigvec = rayleigh_quotient_iteration( eigval, eigvec = rayleigh_quotient_iteration(
linear_operator, case.linear_operator,
numpy.array([0.0, 1.0, 1e-6, 0.0], dtype=complex), case.guess_default,
iterations=8, iterations=8,
) )
residual = norm(operator @ eigvec - eigval * eigvec) residual = norm(case.operator @ eigvec - eigval * eigvec)
assert abs(eigval - 3.0) < 1e-12 assert abs(eigval - 3.0) < 1e-12
assert residual < 1e-12 assert residual < 1e-12
def test_signed_eigensolve_linear_operator_fallback_returns_dominant_positive_mode() -> None: def test_signed_eigensolve_linear_operator_fallback_returns_dominant_positive_mode() -> None:
operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr() case = diagonal_eigen_case()
linear_operator = spalg.LinearOperator(
shape=operator.shape,
dtype=complex,
matvec=lambda vv: operator @ vv,
)
eigvals, eigvecs = signed_eigensolve(linear_operator, how_many=1) eigvals, eigvecs = signed_eigensolve(case.linear_operator, how_many=1)
assert eigvals.shape == (1,) assert eigvals.shape == (1,)
assert eigvecs.shape == (4, 1) assert eigvecs.shape == (4, 1)
assert_allclose(eigvals[0], 5.0, atol=1e-12, rtol=1e-12) assert_close(eigvals[0], 5.0, atol=1e-12, rtol=1e-12)
assert abs(eigvecs[0, 0]) > 0.99 assert abs(eigvecs[0, 0]) > 0.99
def test_power_iteration_finds_dominant_mode() -> None:
case = diagonal_eigen_case()
eigval, eigvec = power_iteration(case.operator, guess_vector=numpy.ones(4, dtype=complex), iterations=20)
assert eigval == pytest.approx(5.0, rel=1e-6)
assert abs(eigvec[0]) > abs(eigvec[1])
def test_rayleigh_quotient_iteration_refines_known_sparse_mode() -> None:
case = diagonal_eigen_case()
def solver(matrix, rhs: numpy.ndarray) -> numpy.ndarray:
return numpy.linalg.lstsq(matrix.toarray(), rhs, rcond=None)[0]
eigval, eigvec = rayleigh_quotient_iteration(
case.operator,
case.guess_sparse,
iterations=8,
solver=solver,
)
residual = numpy.linalg.norm(case.operator @ eigvec - eigval * eigvec)
assert eigval == pytest.approx(3.0, rel=1e-6)
assert residual < 1e-8
def test_signed_eigensolve_returns_largest_positive_modes() -> None:
case = diagonal_eigen_case()
eigvals, eigvecs = signed_eigensolve(case.operator, how_many=2)
assert_close(eigvals, [3.0, 5.0], atol=1e-6)
assert eigvecs.shape == (4, 2)

63
meanas/test/test_eme_numerics.py
View file

@ -1,20 +1,21 @@
import numpy import numpy
from numpy.testing import assert_allclose
from scipy import sparse from scipy import sparse
from ..fdmath import vec from ..fdmath import vec
from ..fdfd import eme from ..fdfd import eme
from ._test_builders import complex_ramp, unit_dxes
from .utils import assert_close
SHAPE = (3, 2, 2) SHAPE = (3, 2, 2)
DXES = [[numpy.ones(2), numpy.ones(2)] for _ in range(2)] DXES = unit_dxes((2, 2))
WAVENUMBERS_L = numpy.array([1.0, 0.8]) WAVENUMBERS_L = numpy.array([1.0, 0.8])
WAVENUMBERS_R = numpy.array([0.9, 0.7]) WAVENUMBERS_R = numpy.array([0.9, 0.7])
def _mode(scale: float) -> tuple[numpy.ndarray, numpy.ndarray]: def _mode(scale: float) -> tuple[numpy.ndarray, numpy.ndarray]:
e_field = (numpy.arange(12).reshape(SHAPE) + 1.0 + scale).astype(complex) e_field = complex_ramp(SHAPE, offset=1.0 + scale)
h_field = (numpy.arange(12).reshape(SHAPE) * 0.2 + 2.0 + 0.05j * scale).astype(complex) h_field = complex_ramp(SHAPE, scale=0.2, offset=2.0, imag_offset=0.05 * scale)
return vec(e_field), vec(h_field) return vec(e_field), vec(h_field)
@ -24,6 +25,29 @@ def _mode_sets() -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], list[tuple[
return left_modes, right_modes return left_modes, right_modes
def _gain_only_tr(*args, **kwargs) -> tuple[numpy.ndarray, numpy.ndarray]:
return numpy.array([[2.0, 0.0], [0.0, 0.5]]), numpy.zeros((2, 2))
def _gain_and_reflection_tr(*args, **kwargs) -> tuple[numpy.ndarray, numpy.ndarray]:
return numpy.array([[2.0, 0.0], [0.0, 0.5]]), numpy.array([[0.0, 1.0], [2.0, 0.0]])
def _nonsymmetric_tr(left_marker: object):
def fake_get_tr(_eh_left, wavenumbers_left, _eh_right, _wavenumbers_right, **kwargs):
if wavenumbers_left is left_marker:
return (
numpy.array([[1.0, 2.0], [0.5, 1.0]]),
numpy.array([[0.0, 1.0], [2.0, 0.0]]),
)
return (
numpy.array([[1.0, -1.0], [0.0, 1.0]]),
numpy.array([[0.0, 0.5], [1.5, 0.0]]),
)
return fake_get_tr
def test_get_tr_returns_finite_bounded_transfer_matrices() -> None: def test_get_tr_returns_finite_bounded_transfer_matrices() -> None:
left_modes, right_modes = _mode_sets() left_modes, right_modes = _mode_sets()
@ -58,7 +82,7 @@ def test_get_abcd_matches_explicit_block_formula() -> None:
assert sparse.issparse(abcd) assert sparse.issparse(abcd)
assert abcd.shape == (4, 4) assert abcd.shape == (4, 4)
assert_allclose(abcd.toarray(), expected) assert_close(abcd.toarray(), expected)
def test_get_s_plain_matches_block_assembly_from_get_tr() -> None: def test_get_s_plain_matches_block_assembly_from_get_tr() -> None:
@ -71,14 +95,11 @@ def test_get_s_plain_matches_block_assembly_from_get_tr() -> None:
assert ss.shape == (4, 4) assert ss.shape == (4, 4)
assert numpy.isfinite(ss).all() assert numpy.isfinite(ss).all()
assert_allclose(ss, expected) assert_close(ss, expected)
def test_get_s_force_nogain_caps_singular_values(monkeypatch) -> None: def test_get_s_force_nogain_caps_singular_values(monkeypatch) -> None:
def fake_get_tr(*args, **kwargs): monkeypatch.setattr(eme, 'get_tr', _gain_only_tr)
return numpy.array([[2.0, 0.0], [0.0, 0.5]]), numpy.zeros((2, 2))
monkeypatch.setattr(eme, 'get_tr', fake_get_tr)
plain_s = eme.get_s(None, None, None, None) plain_s = eme.get_s(None, None, None, None)
clipped_s = eme.get_s(None, None, None, None, force_nogain=True) clipped_s = eme.get_s(None, None, None, None, force_nogain=True)
@ -95,31 +116,17 @@ def test_get_s_force_reciprocal_symmetrizes_output(monkeypatch) -> None:
left = object() left = object()
right = object() right = object()
def fake_get_tr(_eh_left, wavenumbers_left, _eh_right, _wavenumbers_right, **kwargs): monkeypatch.setattr(eme, 'get_tr', _nonsymmetric_tr(left))
if wavenumbers_left is left:
return (
numpy.array([[1.0, 2.0], [0.5, 1.0]]),
numpy.array([[0.0, 1.0], [2.0, 0.0]]),
)
return (
numpy.array([[1.0, -1.0], [0.0, 1.0]]),
numpy.array([[0.0, 0.5], [1.5, 0.0]]),
)
monkeypatch.setattr(eme, 'get_tr', fake_get_tr)
ss = eme.get_s(None, left, None, right, force_reciprocal=True) ss = eme.get_s(None, left, None, right, force_reciprocal=True)
assert_allclose(ss, ss.T) assert_close(ss, ss.T)
def test_get_s_force_nogain_and_reciprocal_returns_finite_output(monkeypatch) -> None: def test_get_s_force_nogain_and_reciprocal_returns_finite_output(monkeypatch) -> None:
def fake_get_tr(*args, **kwargs): monkeypatch.setattr(eme, 'get_tr', _gain_and_reflection_tr)
return numpy.array([[2.0, 0.0], [0.0, 0.5]]), numpy.array([[0.0, 1.0], [2.0, 0.0]])
monkeypatch.setattr(eme, 'get_tr', fake_get_tr)
ss = eme.get_s(None, None, None, None, force_nogain=True, force_reciprocal=True) ss = eme.get_s(None, None, None, None, force_nogain=True, force_reciprocal=True)
assert ss.shape == (4, 4) assert ss.shape == (4, 4)
assert numpy.isfinite(ss).all() assert numpy.isfinite(ss).all()
assert_allclose(ss, ss.T) assert_close(ss, ss.T)
assert (numpy.linalg.svd(ss, compute_uv=False) <= 1.0 + 1e-12).all() assert (numpy.linalg.svd(ss, compute_uv=False) <= 1.0 + 1e-12).all()

122
meanas/test/test_fdfd_algebra_helpers.py
View file

@ -1,52 +1,25 @@
import numpy import numpy
from numpy.testing import assert_allclose
from ..fdmath import vec, unvec from ..fdmath import vec, unvec
from ..fdmath import functional as fd_functional from ..fdmath import functional as fd_functional
from ..fdfd import operators, scpml from ..fdfd import operators, scpml
from ._fdfd_case import (
BOUNDARY_DXES,
OMEGA = 1 / 1500 BOUNDARY_EPSILON,
SHAPE = (2, 3, 2) BOUNDARY_FIELD,
DXES = [ BOUNDARY_SHAPE,
[numpy.array([1.0, 1.5]), numpy.array([0.75, 1.25, 1.5]), numpy.array([1.2, 0.8])], DXES,
[numpy.array([0.9, 1.4]), numpy.array([0.8, 1.1, 1.4]), numpy.array([1.0, 0.7])], EPSILON,
] E_FIELD,
MU,
EPSILON = numpy.stack([ H_FIELD,
numpy.linspace(1.0, 2.2, numpy.prod(SHAPE)).reshape(SHAPE), OMEGA,
numpy.linspace(1.1, 2.3, numpy.prod(SHAPE)).reshape(SHAPE), PEC,
numpy.linspace(1.2, 2.4, numpy.prod(SHAPE)).reshape(SHAPE), PMC,
]) SHAPE,
MU = numpy.stack([ apply_fdfd_matrix,
numpy.linspace(2.0, 3.2, numpy.prod(SHAPE)).reshape(SHAPE), )
numpy.linspace(2.1, 3.3, numpy.prod(SHAPE)).reshape(SHAPE), from .utils import assert_close, assert_fields_close
numpy.linspace(2.2, 3.4, numpy.prod(SHAPE)).reshape(SHAPE),
])
H_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) * 0.25 - 0.75j).astype(complex)
E_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) + 0.5j).astype(complex)
PEC = numpy.zeros((3, *SHAPE), dtype=float)
PEC[1, 0, 1, 0] = 1.0
PMC = numpy.zeros((3, *SHAPE), dtype=float)
PMC[2, 1, 2, 1] = 1.0
BOUNDARY_SHAPE = (3, 4, 3)
BOUNDARY_DXES = [
[numpy.array([1.0, 1.5, 0.8]), numpy.array([0.75, 1.25, 1.5, 0.9]), numpy.array([1.2, 0.8, 1.1])],
[numpy.array([0.9, 1.4, 1.0]), numpy.array([0.8, 1.1, 1.4, 1.0]), numpy.array([1.0, 0.7, 1.3])],
]
BOUNDARY_EPSILON = numpy.stack([
numpy.linspace(1.0, 2.2, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
numpy.linspace(1.1, 2.3, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
numpy.linspace(1.2, 2.4, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
])
BOUNDARY_FIELD = (numpy.arange(3 * numpy.prod(BOUNDARY_SHAPE)).reshape((3, *BOUNDARY_SHAPE)) + 0.5j).astype(complex)
def _apply_matrix(op: operators.sparse.spmatrix, field: numpy.ndarray, shape: tuple[int, ...]) -> numpy.ndarray:
return unvec(op @ vec(field), shape)
def _dense_e_full(mu: numpy.ndarray | None) -> numpy.ndarray: def _dense_e_full(mu: numpy.ndarray | None) -> numpy.ndarray:
@ -84,36 +57,33 @@ def _normalized_distance(u: numpy.ndarray, size: int, thickness: int) -> numpy.n
def test_h_full_matches_dense_reference_with_and_without_mu() -> None: def test_h_full_matches_dense_reference_with_and_without_mu() -> None:
for mu in (None, MU): for mu in (None, MU):
matrix_result = _apply_matrix( matrix_result = apply_fdfd_matrix(
operators.h_full(OMEGA, DXES, vec(EPSILON), None if mu is None else vec(mu), vec(PEC), vec(PMC)), operators.h_full(OMEGA, DXES, vec(EPSILON), None if mu is None else vec(mu), vec(PEC), vec(PMC)),
H_FIELD, H_FIELD,
SHAPE,
) )
dense_result = _dense_h_full(mu) dense_result = _dense_h_full(mu)
assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10) assert_fields_close(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
def test_e_full_matches_dense_reference_with_masks() -> None: def test_e_full_matches_dense_reference_with_masks() -> None:
for mu in (None, MU): for mu in (None, MU):
matrix_result = _apply_matrix( matrix_result = apply_fdfd_matrix(
operators.e_full(OMEGA, DXES, vec(EPSILON), None if mu is None else vec(mu), vec(PEC), vec(PMC)), operators.e_full(OMEGA, DXES, vec(EPSILON), None if mu is None else vec(mu), vec(PEC), vec(PMC)),
E_FIELD, E_FIELD,
SHAPE,
) )
dense_result = _dense_e_full(mu) dense_result = _dense_e_full(mu)
assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10) assert_fields_close(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
def test_h_full_without_masks_matches_dense_reference() -> None: def test_h_full_without_masks_matches_dense_reference() -> None:
ce = fd_functional.curl_forward(DXES[0]) ce = fd_functional.curl_forward(DXES[0])
ch = fd_functional.curl_back(DXES[1]) ch = fd_functional.curl_back(DXES[1])
dense_result = ce(ch(H_FIELD) / EPSILON) - OMEGA**2 * MU * H_FIELD dense_result = ce(ch(H_FIELD) / EPSILON) - OMEGA**2 * MU * H_FIELD
matrix_result = _apply_matrix( matrix_result = apply_fdfd_matrix(
operators.h_full(OMEGA, DXES, vec(EPSILON), vec(MU)), operators.h_full(OMEGA, DXES, vec(EPSILON), vec(MU)),
H_FIELD, H_FIELD,
SHAPE,
) )
assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10) assert_fields_close(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
def test_eh_full_matches_manual_block_operator_with_masks() -> None: def test_eh_full_matches_manual_block_operator_with_masks() -> None:
@ -130,23 +100,23 @@ def test_eh_full_matches_manual_block_operator_with_masks() -> None:
dense_e = pe * ch(pm * H_FIELD) - pe * (1j * OMEGA * EPSILON * (pe * E_FIELD)) dense_e = pe * ch(pm * H_FIELD) - pe * (1j * OMEGA * EPSILON * (pe * E_FIELD))
dense_h = pm * ce(pe * E_FIELD) + pm * (1j * OMEGA * MU * (pm * H_FIELD)) dense_h = pm * ce(pe * E_FIELD) + pm * (1j * OMEGA * MU * (pm * H_FIELD))
assert_allclose(matrix_e, dense_e, atol=1e-10, rtol=1e-10) assert_fields_close(matrix_e, dense_e, atol=1e-10, rtol=1e-10)
assert_allclose(matrix_h, dense_h, atol=1e-10, rtol=1e-10) assert_fields_close(matrix_h, dense_h, atol=1e-10, rtol=1e-10)
def test_e2h_pmc_mask_matches_masked_unmasked_result() -> None: def test_e2h_pmc_mask_matches_masked_unmasked_result() -> None:
pmc_complement = numpy.where(PMC, 0.0, 1.0) pmc_complement = numpy.where(PMC, 0.0, 1.0)
unmasked = _apply_matrix(operators.e2h(OMEGA, DXES, vec(MU)), E_FIELD, SHAPE) unmasked = apply_fdfd_matrix(operators.e2h(OMEGA, DXES, vec(MU)), E_FIELD)
masked = _apply_matrix(operators.e2h(OMEGA, DXES, vec(MU), vec(PMC)), E_FIELD, SHAPE) masked = apply_fdfd_matrix(operators.e2h(OMEGA, DXES, vec(MU), vec(PMC)), E_FIELD)
assert_allclose(masked, pmc_complement * unmasked, atol=1e-10, rtol=1e-10) assert_fields_close(masked, pmc_complement * unmasked, atol=1e-10, rtol=1e-10)
def test_poynting_h_cross_matches_negative_e_cross_relation() -> None: def test_poynting_h_cross_matches_negative_e_cross_relation() -> None:
h_cross_e = _apply_matrix(operators.poynting_h_cross(vec(H_FIELD), DXES), E_FIELD, SHAPE) h_cross_e = apply_fdfd_matrix(operators.poynting_h_cross(vec(H_FIELD), DXES), E_FIELD)
e_cross_h = _apply_matrix(operators.poynting_e_cross(vec(E_FIELD), DXES), H_FIELD, SHAPE) e_cross_h = apply_fdfd_matrix(operators.poynting_e_cross(vec(E_FIELD), DXES), H_FIELD)
assert_allclose(h_cross_e, -e_cross_h, atol=1e-10, rtol=1e-10) assert_fields_close(h_cross_e, -e_cross_h, atol=1e-10, rtol=1e-10)
def test_e_boundary_source_interior_mask_is_independent_of_periodic_edges() -> None: def test_e_boundary_source_interior_mask_is_independent_of_periodic_edges() -> None:
@ -156,7 +126,7 @@ def test_e_boundary_source_interior_mask_is_independent_of_periodic_edges() -> N
periodic = operators.e_boundary_source(vec(mask), OMEGA, BOUNDARY_DXES, vec(BOUNDARY_EPSILON), periodic_mask_edges=True) periodic = operators.e_boundary_source(vec(mask), OMEGA, BOUNDARY_DXES, vec(BOUNDARY_EPSILON), periodic_mask_edges=True)
mirrored = operators.e_boundary_source(vec(mask), OMEGA, BOUNDARY_DXES, vec(BOUNDARY_EPSILON), periodic_mask_edges=False) mirrored = operators.e_boundary_source(vec(mask), OMEGA, BOUNDARY_DXES, vec(BOUNDARY_EPSILON), periodic_mask_edges=False)
assert_allclose(periodic.toarray(), mirrored.toarray()) assert_close(periodic.toarray(), mirrored.toarray())
def test_e_boundary_source_periodic_edges_add_opposite_face_response() -> None: def test_e_boundary_source_periodic_edges_add_opposite_face_response() -> None:
@ -168,7 +138,7 @@ def test_e_boundary_source_periodic_edges_add_opposite_face_response() -> None:
diff = unvec((periodic - mirrored) @ vec(BOUNDARY_FIELD), BOUNDARY_SHAPE) diff = unvec((periodic - mirrored) @ vec(BOUNDARY_FIELD), BOUNDARY_SHAPE)
assert numpy.isfinite(diff).all() assert numpy.isfinite(diff).all()
assert_allclose(diff[:, 1:-1, :, :], 0.0) assert_close(diff[:, 1:-1, :, :], 0.0)
assert numpy.linalg.norm(diff[:, -1, :, :]) > 0 assert numpy.linalg.norm(diff[:, -1, :, :]) > 0
@ -179,7 +149,7 @@ def test_prepare_s_function_matches_closed_form_polynomial() -> None:
s_function = scpml.prepare_s_function(ln_R=ln_r, m=order) s_function = scpml.prepare_s_function(ln_R=ln_r, m=order)
expected = (order + 1) * ln_r / 2 * distances**order expected = (order + 1) * ln_r / 2 * distances**order
assert_allclose(s_function(distances), expected) assert_close(s_function(distances), expected)
def test_uniform_grid_scpml_matches_expected_stretch_profile() -> None: def test_uniform_grid_scpml_matches_expected_stretch_profile() -> None:
@ -191,11 +161,11 @@ def test_uniform_grid_scpml_matches_expected_stretch_profile() -> None:
grid = numpy.arange(size, dtype=float) grid = numpy.arange(size, dtype=float)
expected_a = 1 + 1j * s_function(_normalized_distance(grid, size, thickness)) / correction expected_a = 1 + 1j * s_function(_normalized_distance(grid, size, thickness)) / correction
expected_b = 1 + 1j * s_function(_normalized_distance(grid + 0.5, size, thickness)) / correction expected_b = 1 + 1j * s_function(_normalized_distance(grid + 0.5, size, thickness)) / correction
assert_allclose(dxes[0][axis], expected_a) assert_close(dxes[0][axis], expected_a)
assert_allclose(dxes[1][axis], expected_b) assert_close(dxes[1][axis], expected_b)
assert_allclose(dxes[0][1], 1.0) assert_close(dxes[0][1], 1.0)
assert_allclose(dxes[1][1], 1.0) assert_close(dxes[1][1], 1.0)
assert numpy.isfinite(dxes[0][0]).all() assert numpy.isfinite(dxes[0][0]).all()
assert numpy.isfinite(dxes[1][0]).all() assert numpy.isfinite(dxes[1][0]).all()
@ -206,7 +176,7 @@ def test_uniform_grid_scpml_default_s_function_matches_explicit_default() -> Non
for implicit_group, explicit_group in zip(implicit, explicit, strict=True): for implicit_group, explicit_group in zip(implicit, explicit, strict=True):
for implicit_axis, explicit_axis in zip(implicit_group, explicit_group, strict=True): for implicit_axis, explicit_axis in zip(implicit_group, explicit_group, strict=True):
assert_allclose(implicit_axis, explicit_axis) assert_close(implicit_axis, explicit_axis)
def test_stretch_with_scpml_only_modifies_requested_front_edge() -> None: def test_stretch_with_scpml_only_modifies_requested_front_edge() -> None:
@ -214,10 +184,10 @@ def test_stretch_with_scpml_only_modifies_requested_front_edge() -> None:
base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)] base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)]
stretched = scpml.stretch_with_scpml(base, axis=0, polarity=1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function) stretched = scpml.stretch_with_scpml(base, axis=0, polarity=1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function)
assert_allclose(stretched[0][0][2:], 1.0) assert_close(stretched[0][0][2:], 1.0)
assert_allclose(stretched[1][0][2:], 1.0) assert_close(stretched[1][0][2:], 1.0)
assert_allclose(stretched[0][0][-2:], 1.0) assert_close(stretched[0][0][-2:], 1.0)
assert_allclose(stretched[1][0][-2:], 1.0) assert_close(stretched[1][0][-2:], 1.0)
assert numpy.linalg.norm(stretched[0][0][:2] - 1.0) > 0 assert numpy.linalg.norm(stretched[0][0][:2] - 1.0) > 0
assert numpy.linalg.norm(stretched[1][0][:2] - 1.0) > 0 assert numpy.linalg.norm(stretched[1][0][:2] - 1.0) > 0
@ -227,8 +197,8 @@ def test_stretch_with_scpml_only_modifies_requested_back_edge() -> None:
base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)] base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)]
stretched = scpml.stretch_with_scpml(base, axis=0, polarity=-1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function) stretched = scpml.stretch_with_scpml(base, axis=0, polarity=-1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function)
assert_allclose(stretched[0][0][:4], 1.0) assert_close(stretched[0][0][:4], 1.0)
assert_allclose(stretched[1][0][:4], 1.0) assert_close(stretched[1][0][:4], 1.0)
assert numpy.linalg.norm(stretched[0][0][-2:] - 1.0) > 0 assert numpy.linalg.norm(stretched[0][0][-2:] - 1.0) > 0
assert numpy.linalg.norm(stretched[1][0][-2:] - 1.0) > 0 assert numpy.linalg.norm(stretched[1][0][-2:] - 1.0) > 0
@ -240,4 +210,4 @@ def test_stretch_with_scpml_thickness_zero_is_noop() -> None:
for grid_group in stretched: for grid_group in stretched:
for axis_grid in grid_group: for axis_grid in grid_group:
assert_allclose(axis_grid, 1.0) assert_close(axis_grid, 1.0)

26
meanas/test/test_fdfd_farfield.py
View file

@ -72,3 +72,29 @@ def test_far_to_nearfield_uses_default_and_scalar_padding_shapes() -> None:
assert default_result['H'][0].shape == (8, 8) assert default_result['H'][0].shape == (8, 8)
assert scalar_result['E'][0].shape == (4, 4) assert scalar_result['E'][0].shape == (4, 4)
assert scalar_result['H'][0].shape == (4, 4) assert scalar_result['H'][0].shape == (4, 4)
def test_farfield_roundtrip_supports_rectangular_arrays() -> None:
e_near = [numpy.zeros((4, 8), dtype=complex), numpy.zeros((4, 8), dtype=complex)]
h_near = [numpy.zeros((4, 8), dtype=complex), numpy.zeros((4, 8), dtype=complex)]
e_near[0][1, 3] = 1.0 + 0.25j
h_near[1][2, 5] = -0.5j
ff = farfield.near_to_farfield(e_near, h_near, dx=0.2, dy=0.3, padded_size=(4, 8))
restored = farfield.far_to_nearfield(
[field.copy() for field in ff['E']],
[field.copy() for field in ff['H']],
ff['dkx'],
ff['dky'],
padded_size=(4, 8),
)
assert isinstance(ff['dkx'], float)
assert isinstance(ff['dky'], float)
assert ff['E'][0].shape == (4, 8)
assert restored['E'][0].shape == (4, 8)
assert restored['H'][0].shape == (4, 8)
assert restored['dx'] == pytest.approx(0.2)
assert restored['dy'] == pytest.approx(0.3)
assert numpy.isfinite(restored['E'][0]).all()
assert numpy.isfinite(restored['H'][0]).all()

53
meanas/test/test_fdfd_functional.py
View file

@ -1,48 +1,21 @@
import numpy import numpy
from numpy.testing import assert_allclose
from ..fdmath import vec, unvec from ..fdmath import unvec, vec
from ..fdfd import functional, operators from ..fdfd import functional, operators
from ._fdfd_case import DXES, EPSILON, E_FIELD, H_FIELD, MU, OMEGA, SHAPE, TF_REGION, apply_fdfd_matrix
from .utils import assert_fields_close
OMEGA = 1 / 1500
SHAPE = (2, 3, 2)
ATOL = 1e-9 ATOL = 1e-9
RTOL = 1e-9 RTOL = 1e-9
DXES = [
[numpy.array([1.0, 1.5]), numpy.array([0.75, 1.25, 1.5]), numpy.array([1.2, 0.8])],
[numpy.array([0.9, 1.4]), numpy.array([0.8, 1.1, 1.4]), numpy.array([1.0, 0.7])],
]
EPSILON = numpy.stack([
numpy.linspace(1.0, 2.2, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(1.1, 2.3, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(1.2, 2.4, numpy.prod(SHAPE)).reshape(SHAPE),
])
MU = numpy.stack([
numpy.linspace(2.0, 3.2, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(2.1, 3.3, numpy.prod(SHAPE)).reshape(SHAPE),
numpy.linspace(2.2, 3.4, numpy.prod(SHAPE)).reshape(SHAPE),
])
E_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) + 0.5j).astype(complex)
H_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) * 0.25 - 0.75j).astype(complex)
TF_REGION = numpy.zeros((3, *SHAPE), dtype=float)
TF_REGION[:, 0, 1, 0] = 1.0
def apply_matrix(op: operators.sparse.spmatrix, field: numpy.ndarray) -> numpy.ndarray:
return unvec(op @ vec(field), SHAPE)
def assert_fields_match(actual: numpy.ndarray, expected: numpy.ndarray) -> None: def assert_fields_match(actual: numpy.ndarray, expected: numpy.ndarray) -> None:
assert_allclose(actual, expected, atol=ATOL, rtol=RTOL) assert_fields_close(actual, expected, atol=ATOL, rtol=RTOL)
def test_e_full_matches_sparse_operator_without_mu() -> None: def test_e_full_matches_sparse_operator_without_mu() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.e_full(OMEGA, DXES, vec(EPSILON)), operators.e_full(OMEGA, DXES, vec(EPSILON)),
E_FIELD, E_FIELD,
) )
@ -52,7 +25,7 @@ def test_e_full_matches_sparse_operator_without_mu() -> None:
def test_e_full_matches_sparse_operator_with_mu() -> None: def test_e_full_matches_sparse_operator_with_mu() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.e_full(OMEGA, DXES, vec(EPSILON), vec(MU)), operators.e_full(OMEGA, DXES, vec(EPSILON), vec(MU)),
E_FIELD, E_FIELD,
) )
@ -80,7 +53,7 @@ def test_eh_full_matches_sparse_operator_without_mu() -> None:
def test_e2h_matches_sparse_operator_with_mu() -> None: def test_e2h_matches_sparse_operator_with_mu() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.e2h(OMEGA, DXES, vec(MU)), operators.e2h(OMEGA, DXES, vec(MU)),
E_FIELD, E_FIELD,
) )
@ -90,7 +63,7 @@ def test_e2h_matches_sparse_operator_with_mu() -> None:
def test_e2h_matches_sparse_operator_without_mu() -> None: def test_e2h_matches_sparse_operator_without_mu() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.e2h(OMEGA, DXES), operators.e2h(OMEGA, DXES),
E_FIELD, E_FIELD,
) )
@ -100,7 +73,7 @@ def test_e2h_matches_sparse_operator_without_mu() -> None:
def test_m2j_matches_sparse_operator_without_mu() -> None: def test_m2j_matches_sparse_operator_without_mu() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.m2j(OMEGA, DXES), operators.m2j(OMEGA, DXES),
H_FIELD, H_FIELD,
) )
@ -110,7 +83,7 @@ def test_m2j_matches_sparse_operator_without_mu() -> None:
def test_m2j_matches_sparse_operator_with_mu() -> None: def test_m2j_matches_sparse_operator_with_mu() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.m2j(OMEGA, DXES, vec(MU)), operators.m2j(OMEGA, DXES, vec(MU)),
H_FIELD, H_FIELD,
) )
@ -120,7 +93,7 @@ def test_m2j_matches_sparse_operator_with_mu() -> None:
def test_e_tfsf_source_matches_sparse_operator_without_mu() -> None: def test_e_tfsf_source_matches_sparse_operator_without_mu() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.e_tfsf_source(vec(TF_REGION), OMEGA, DXES, vec(EPSILON)), operators.e_tfsf_source(vec(TF_REGION), OMEGA, DXES, vec(EPSILON)),
E_FIELD, E_FIELD,
) )
@ -130,7 +103,7 @@ def test_e_tfsf_source_matches_sparse_operator_without_mu() -> None:
def test_e_tfsf_source_matches_sparse_operator_with_mu() -> None: def test_e_tfsf_source_matches_sparse_operator_with_mu() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.e_tfsf_source(vec(TF_REGION), OMEGA, DXES, vec(EPSILON), vec(MU)), operators.e_tfsf_source(vec(TF_REGION), OMEGA, DXES, vec(EPSILON), vec(MU)),
E_FIELD, E_FIELD,
) )
@ -140,7 +113,7 @@ def test_e_tfsf_source_matches_sparse_operator_with_mu() -> None:
def test_poynting_e_cross_h_matches_sparse_operator() -> None: def test_poynting_e_cross_h_matches_sparse_operator() -> None:
matrix_result = apply_matrix( matrix_result = apply_fdfd_matrix(
operators.poynting_e_cross(vec(E_FIELD), DXES), operators.poynting_e_cross(vec(E_FIELD), DXES),
H_FIELD, H_FIELD,
) )

106
meanas/test/test_fdfd_solvers.py
View file

@ -1,140 +1,126 @@
import numpy import numpy
from numpy.testing import assert_allclose
from scipy import sparse
from ..fdfd import solvers from ..fdfd import solvers
from ._solver_cases import solver_plumbing_case
from .utils import assert_close
def test_scipy_qmr_wraps_user_callback_without_recursion(monkeypatch) -> None: def test_scipy_qmr_wraps_user_callback_without_recursion(monkeypatch) -> None:
seen: list[tuple[float, ...]] = [] seen: list[tuple[float, ...]] = []
def fake_qmr(a: sparse.spmatrix, b: numpy.ndarray, **kwargs): def fake_qmr(a, b: numpy.ndarray, **kwargs):
kwargs['callback'](numpy.array([1.0, 2.0])) kwargs['callback'](numpy.array([1.0, 2.0]))
return numpy.array([3.0, 4.0]), 0 return numpy.array([3.0, 4.0]), 0
monkeypatch.setattr(solvers.scipy.sparse.linalg, 'qmr', fake_qmr) monkeypatch.setattr(solvers.scipy.sparse.linalg, 'qmr', fake_qmr)
result = solvers._scipy_qmr( result = solvers._scipy_qmr(
sparse.eye(2).tocsr(), solver_plumbing_case().a0,
numpy.array([1.0, 0.0]), numpy.array([1.0, 0.0]),
callback=lambda xk: seen.append(tuple(xk)), callback=lambda xk: seen.append(tuple(xk)),
) )
assert_allclose(result, [3.0, 4.0]) assert_close(result, [3.0, 4.0])
assert seen == [(1.0, 2.0)] assert seen == [(1.0, 2.0)]
def test_scipy_qmr_installs_logging_callback_when_missing(monkeypatch) -> None: def test_scipy_qmr_installs_logging_callback_when_missing(monkeypatch) -> None:
callback_seen: list[numpy.ndarray] = [] callback_seen: list[numpy.ndarray] = []
def fake_qmr(a: sparse.spmatrix, b: numpy.ndarray, **kwargs): def fake_qmr(a, b: numpy.ndarray, **kwargs):
callback = kwargs['callback'] callback = kwargs['callback']
callback(numpy.array([5.0, 6.0])) callback(numpy.array([5.0, 6.0]))
callback_seen.append(b.copy()) callback_seen.append(b.copy())
return numpy.array([7.0, 8.0]), 0 return numpy.array([7.0, 8.0]), 0
monkeypatch.setattr(solvers.scipy.sparse.linalg, 'qmr', fake_qmr) monkeypatch.setattr(solvers.scipy.sparse.linalg, 'qmr', fake_qmr)
result = solvers._scipy_qmr(sparse.eye(2).tocsr(), numpy.array([1.0, 0.0])) result = solvers._scipy_qmr(solver_plumbing_case().a0, numpy.array([1.0, 0.0]))
assert_allclose(result, [7.0, 8.0]) assert_close(result, [7.0, 8.0])
assert len(callback_seen) == 1 assert len(callback_seen) == 1
def test_generic_forward_preconditions_system_and_guess(monkeypatch) -> None: def test_generic_forward_preconditions_system_and_guess(monkeypatch) -> None:
omega = 2.0 case = solver_plumbing_case()
a0 = sparse.csr_matrix(numpy.array([[1.0 + 2.0j, 2.0], [3.0 - 1.0j, 4.0]])) captured: dict[str, object] = {}
pl = sparse.csr_matrix(numpy.array([[2.0, 0.0], [0.0, 3.0j]]))
pr = sparse.csr_matrix(numpy.array([[0.5, 0.0], [0.0, -2.0j]]))
j = numpy.array([1.0 + 0.5j, -2.0])
guess = numpy.array([0.25 - 0.75j, 1.5 + 0.5j])
solver_result = numpy.array([3.0 - 1.0j, -4.0 + 2.0j])
captured: dict[str, numpy.ndarray | sparse.spmatrix] = {}
monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: a0) monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: case.a0)
monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (pl, pr)) monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (case.pl, case.pr))
def fake_solver(a: sparse.spmatrix, b: numpy.ndarray, **kwargs): def fake_solver(a, b: numpy.ndarray, **kwargs):
captured['a'] = a captured['a'] = a
captured['b'] = b captured['b'] = b
captured['x0'] = kwargs['x0'] captured['x0'] = kwargs['x0']
captured['atol'] = kwargs['atol'] captured['atol'] = kwargs['atol']
return solver_result return case.solver_result
result = solvers.generic( result = solvers.generic(
omega=omega, omega=case.omega,
dxes=[[numpy.ones(1) for _ in range(3)] for _ in range(2)], dxes=case.dxes,
J=j, J=case.j,
epsilon=numpy.ones(2), epsilon=case.epsilon,
matrix_solver=fake_solver, matrix_solver=fake_solver,
matrix_solver_opts={'atol': 1e-12}, matrix_solver_opts={'atol': 1e-12},
E_guess=guess, E_guess=case.guess,
) )
assert_allclose(captured['a'].toarray(), (pl @ a0 @ pr).toarray()) assert_close(captured['a'].toarray(), (case.pl @ case.a0 @ case.pr).toarray())
assert_allclose(captured['b'], pl @ (-1j * omega * j)) assert_close(captured['b'], case.pl @ (-1j * case.omega * case.j))
assert_allclose(captured['x0'], pl @ guess) assert_close(captured['x0'], case.pl @ case.guess)
assert captured['atol'] == 1e-12 assert captured['atol'] == 1e-12
assert_allclose(result, pr @ solver_result) assert_close(result, case.pr @ case.solver_result)
def test_generic_adjoint_preconditions_system_and_guess(monkeypatch) -> None: def test_generic_adjoint_preconditions_system_and_guess(monkeypatch) -> None:
omega = 2.0 case = solver_plumbing_case()
a0 = sparse.csr_matrix(numpy.array([[1.0 + 2.0j, 2.0], [3.0 - 1.0j, 4.0]])) captured: dict[str, object] = {}
pl = sparse.csr_matrix(numpy.array([[2.0, 0.0], [0.0, 3.0j]]))
pr = sparse.csr_matrix(numpy.array([[0.5, 0.0], [0.0, -2.0j]]))
j = numpy.array([1.0 + 0.5j, -2.0])
guess = numpy.array([0.25 - 0.75j, 1.5 + 0.5j])
solver_result = numpy.array([3.0 - 1.0j, -4.0 + 2.0j])
captured: dict[str, numpy.ndarray | sparse.spmatrix] = {}
monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: a0) monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: case.a0)
monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (pl, pr)) monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (case.pl, case.pr))
def fake_solver(a: sparse.spmatrix, b: numpy.ndarray, **kwargs): def fake_solver(a, b: numpy.ndarray, **kwargs):
captured['a'] = a captured['a'] = a
captured['b'] = b captured['b'] = b
captured['x0'] = kwargs['x0'] captured['x0'] = kwargs['x0']
captured['rtol'] = kwargs['rtol'] captured['rtol'] = kwargs['rtol']
return solver_result return case.solver_result
result = solvers.generic( result = solvers.generic(
omega=omega, omega=case.omega,
dxes=[[numpy.ones(1) for _ in range(3)] for _ in range(2)], dxes=case.dxes,
J=j, J=case.j,
epsilon=numpy.ones(2), epsilon=case.epsilon,
matrix_solver=fake_solver, matrix_solver=fake_solver,
matrix_solver_opts={'rtol': 1e-9}, matrix_solver_opts={'rtol': 1e-9},
E_guess=guess, E_guess=case.guess,
adjoint=True, adjoint=True,
) )
expected_matrix = (pl @ a0 @ pr).T.conjugate() expected_matrix = (case.pl @ case.a0 @ case.pr).T.conjugate()
assert_allclose(captured['a'].toarray(), expected_matrix.toarray()) assert_close(captured['a'].toarray(), expected_matrix.toarray())
assert_allclose(captured['b'], pr.T.conjugate() @ (-1j * omega * j)) assert_close(captured['b'], case.pr.T.conjugate() @ (-1j * case.omega * case.j))
assert_allclose(captured['x0'], pr.T.conjugate() @ guess) assert_close(captured['x0'], case.pr.T.conjugate() @ case.guess)
assert captured['rtol'] == 1e-9 assert captured['rtol'] == 1e-9
assert_allclose(result, pl.T.conjugate() @ solver_result) assert_close(result, case.pl.T.conjugate() @ case.solver_result)
def test_generic_without_guess_does_not_inject_x0(monkeypatch) -> None: def test_generic_without_guess_does_not_inject_x0(monkeypatch) -> None:
a0 = sparse.eye(2).tocsr() case = solver_plumbing_case()
pl = sparse.eye(2).tocsr()
pr = sparse.eye(2).tocsr()
captured: dict[str, object] = {} captured: dict[str, object] = {}
monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: a0) monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: case.a0)
monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (pl, pr)) monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (case.pl, case.pr))
def fake_solver(a: sparse.spmatrix, b: numpy.ndarray, **kwargs): def fake_solver(a, b: numpy.ndarray, **kwargs):
captured['kwargs'] = kwargs captured['kwargs'] = kwargs
return numpy.array([1.0, -1.0]) return numpy.array([1.0, -1.0])
result = solvers.generic( result = solvers.generic(
omega=1.0, omega=1.0,
dxes=[[numpy.ones(1) for _ in range(3)] for _ in range(2)], dxes=case.dxes,
J=numpy.array([2.0, 3.0]), J=numpy.array([2.0, 3.0]),
epsilon=numpy.ones(2), epsilon=case.epsilon,
matrix_solver=fake_solver, matrix_solver=fake_solver,
) )
assert 'x0' not in captured['kwargs'] assert 'x0' not in captured['kwargs']
assert_allclose(result, [1.0, -1.0]) assert_close(result, case.pr @ numpy.array([1.0, -1.0]))

18
meanas/test/test_fdmath_functional.py
View file

@ -1,9 +1,9 @@
import numpy import numpy
from numpy.testing import assert_allclose
from ..fdmath import functional as fd_functional from ..fdmath import functional as fd_functional
from ..fdmath import operators as fd_operators from ..fdmath import operators as fd_operators
from ..fdmath import vec, unvec from ..fdmath import vec, unvec
from .utils import assert_close, assert_fields_close
SHAPE = (2, 3, 2) SHAPE = (2, 3, 2)
@ -20,13 +20,13 @@ VECTOR_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) + 0.25j
def test_deriv_forward_without_dx_matches_numpy_roll() -> None: def test_deriv_forward_without_dx_matches_numpy_roll() -> None:
for axis, deriv in enumerate(fd_functional.deriv_forward()): for axis, deriv in enumerate(fd_functional.deriv_forward()):
expected = numpy.roll(SCALAR_FIELD, -1, axis=axis) - SCALAR_FIELD expected = numpy.roll(SCALAR_FIELD, -1, axis=axis) - SCALAR_FIELD
assert_allclose(deriv(SCALAR_FIELD), expected) assert_close(deriv(SCALAR_FIELD), expected)
def test_deriv_back_without_dx_matches_numpy_roll() -> None: def test_deriv_back_without_dx_matches_numpy_roll() -> None:
for axis, deriv in enumerate(fd_functional.deriv_back()): for axis, deriv in enumerate(fd_functional.deriv_back()):
expected = SCALAR_FIELD - numpy.roll(SCALAR_FIELD, 1, axis=axis) expected = SCALAR_FIELD - numpy.roll(SCALAR_FIELD, 1, axis=axis)
assert_allclose(deriv(SCALAR_FIELD), expected) assert_close(deriv(SCALAR_FIELD), expected)
def test_curl_parts_sum_to_full_curl() -> None: def test_curl_parts_sum_to_full_curl() -> None:
@ -36,18 +36,18 @@ def test_curl_parts_sum_to_full_curl() -> None:
back_parts = fd_functional.curl_back_parts(DX_H)(VECTOR_FIELD) back_parts = fd_functional.curl_back_parts(DX_H)(VECTOR_FIELD)
for axis in range(3): for axis in range(3):
assert_allclose(forward_parts[axis][0] + forward_parts[axis][1], curl_forward[axis]) assert_close(forward_parts[axis][0] + forward_parts[axis][1], curl_forward[axis])
assert_allclose(back_parts[axis][0] + back_parts[axis][1], curl_back[axis]) assert_close(back_parts[axis][0] + back_parts[axis][1], curl_back[axis])
def test_derivatives_match_sparse_operators_on_nonuniform_grid() -> None: def test_derivatives_match_sparse_operators_on_nonuniform_grid() -> None:
for axis, deriv in enumerate(fd_functional.deriv_forward(DX_E)): for axis, deriv in enumerate(fd_functional.deriv_forward(DX_E)):
matrix_result = (fd_operators.deriv_forward(DX_E)[axis] @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C') matrix_result = (fd_operators.deriv_forward(DX_E)[axis] @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C')
assert_allclose(deriv(SCALAR_FIELD), matrix_result, atol=1e-12, rtol=1e-12) assert_close(deriv(SCALAR_FIELD), matrix_result, atol=1e-12, rtol=1e-12)
for axis, deriv in enumerate(fd_functional.deriv_back(DX_H)): for axis, deriv in enumerate(fd_functional.deriv_back(DX_H)):
matrix_result = (fd_operators.deriv_back(DX_H)[axis] @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C') matrix_result = (fd_operators.deriv_back(DX_H)[axis] @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C')
assert_allclose(deriv(SCALAR_FIELD), matrix_result, atol=1e-12, rtol=1e-12) assert_close(deriv(SCALAR_FIELD), matrix_result, atol=1e-12, rtol=1e-12)
def test_curls_match_sparse_operators_on_nonuniform_grid() -> None: def test_curls_match_sparse_operators_on_nonuniform_grid() -> None:
@ -56,5 +56,5 @@ def test_curls_match_sparse_operators_on_nonuniform_grid() -> None:
matrix_forward = unvec(fd_operators.curl_forward(DX_E) @ vec(VECTOR_FIELD), SHAPE) matrix_forward = unvec(fd_operators.curl_forward(DX_E) @ vec(VECTOR_FIELD), SHAPE)
matrix_back = unvec(fd_operators.curl_back(DX_H) @ vec(VECTOR_FIELD), SHAPE) matrix_back = unvec(fd_operators.curl_back(DX_H) @ vec(VECTOR_FIELD), SHAPE)
assert_allclose(curl_forward, matrix_forward, atol=1e-12, rtol=1e-12) assert_fields_close(curl_forward, matrix_forward, atol=1e-12, rtol=1e-12)
assert_allclose(curl_back, matrix_back, atol=1e-12, rtol=1e-12) assert_fields_close(curl_back, matrix_back, atol=1e-12, rtol=1e-12)

19
meanas/test/test_fdmath_operators.py
View file

@ -1,14 +1,15 @@
import numpy import numpy
import pytest import pytest
from numpy.testing import assert_allclose, assert_array_equal
from ..fdmath import operators, unvec, vec from ..fdmath import operators, unvec, vec
from ._test_builders import real_ramp
from .utils import assert_close
SHAPE = (2, 3, 2) SHAPE = (2, 3, 2)
SCALAR_FIELD = numpy.arange(numpy.prod(SHAPE), dtype=float).reshape(SHAPE, order='C') SCALAR_FIELD = real_ramp(SHAPE)
VECTOR_LEFT = (numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') + 0.5) VECTOR_LEFT = real_ramp((3, *SHAPE), offset=0.5)
VECTOR_RIGHT = (2.0 + numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') / 3.0) VECTOR_RIGHT = real_ramp((3, *SHAPE), scale=1 / 3, offset=2.0)
def _apply_scalar_matrix(op: operators.sparse.spmatrix) -> numpy.ndarray: def _apply_scalar_matrix(op: operators.sparse.spmatrix) -> numpy.ndarray:
@ -26,7 +27,7 @@ def _mirrored_indices(size: int, shift_distance: int) -> numpy.ndarray:
def test_shift_circ_matches_numpy_roll(axis: int, shift_distance: int) -> None: def test_shift_circ_matches_numpy_roll(axis: int, shift_distance: int) -> None:
matrix_result = _apply_scalar_matrix(operators.shift_circ(axis, SHAPE, shift_distance)) matrix_result = _apply_scalar_matrix(operators.shift_circ(axis, SHAPE, shift_distance))
expected = numpy.roll(SCALAR_FIELD, -shift_distance, axis=axis) expected = numpy.roll(SCALAR_FIELD, -shift_distance, axis=axis)
assert_array_equal(matrix_result, expected) assert_close(matrix_result, expected)
@pytest.mark.parametrize(('axis', 'shift_distance'), [(0, 1), (1, -1), (2, 1)]) @pytest.mark.parametrize(('axis', 'shift_distance'), [(0, 1), (1, -1), (2, 1)])
@ -35,7 +36,7 @@ def test_shift_with_mirror_matches_explicit_mirrored_indices(axis: int, shift_di
indices = [numpy.arange(length) for length in SHAPE] indices = [numpy.arange(length) for length in SHAPE]
indices[axis] = _mirrored_indices(SHAPE[axis], shift_distance) indices[axis] = _mirrored_indices(SHAPE[axis], shift_distance)
expected = SCALAR_FIELD[numpy.ix_(*indices)] expected = SCALAR_FIELD[numpy.ix_(*indices)]
assert_array_equal(matrix_result, expected) assert_close(matrix_result, expected)
@pytest.mark.parametrize( @pytest.mark.parametrize(
@ -69,19 +70,19 @@ def test_vec_cross_matches_pointwise_cross_product() -> None:
expected[0] = VECTOR_LEFT[1] * VECTOR_RIGHT[2] - VECTOR_LEFT[2] * VECTOR_RIGHT[1] expected[0] = VECTOR_LEFT[1] * VECTOR_RIGHT[2] - VECTOR_LEFT[2] * VECTOR_RIGHT[1]
expected[1] = VECTOR_LEFT[2] * VECTOR_RIGHT[0] - VECTOR_LEFT[0] * VECTOR_RIGHT[2] expected[1] = VECTOR_LEFT[2] * VECTOR_RIGHT[0] - VECTOR_LEFT[0] * VECTOR_RIGHT[2]
expected[2] = VECTOR_LEFT[0] * VECTOR_RIGHT[1] - VECTOR_LEFT[1] * VECTOR_RIGHT[0] expected[2] = VECTOR_LEFT[0] * VECTOR_RIGHT[1] - VECTOR_LEFT[1] * VECTOR_RIGHT[0]
assert_allclose(matrix_result, expected) assert_close(matrix_result, expected)
def test_avg_forward_matches_half_sum_with_forward_neighbor() -> None: def test_avg_forward_matches_half_sum_with_forward_neighbor() -> None:
matrix_result = _apply_scalar_matrix(operators.avg_forward(1, SHAPE)) matrix_result = _apply_scalar_matrix(operators.avg_forward(1, SHAPE))
expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, -1, axis=1)) expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, -1, axis=1))
assert_allclose(matrix_result, expected) assert_close(matrix_result, expected)
def test_avg_back_matches_half_sum_with_backward_neighbor() -> None: def test_avg_back_matches_half_sum_with_backward_neighbor() -> None:
matrix_result = _apply_scalar_matrix(operators.avg_back(1, SHAPE)) matrix_result = _apply_scalar_matrix(operators.avg_back(1, SHAPE))
expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, 1, axis=1)) expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, 1, axis=1))
assert_allclose(matrix_result, expected) assert_close(matrix_result, expected)
def test_avg_forward_rejects_invalid_shape() -> None: def test_avg_forward_rejects_invalid_shape() -> None:

15
meanas/test/test_fdmath_vectorization.py
View file

@ -1,12 +1,13 @@
import numpy import numpy
from numpy.testing import assert_allclose, assert_array_equal
from ..fdmath import unvec, vec from ..fdmath import unvec, vec
from ._test_builders import complex_ramp, real_ramp
from .utils import assert_close
SHAPE = (2, 3, 2) SHAPE = (2, 3, 2)
FIELD = numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') FIELD = real_ramp((3, *SHAPE))
COMPLEX_FIELD = (FIELD + 0.5j * FIELD).astype(complex) COMPLEX_FIELD = complex_ramp((3, *SHAPE), imag_scale=0.5)
def test_vec_and_unvec_return_none_for_none_input() -> None: def test_vec_and_unvec_return_none_for_none_input() -> None:
@ -20,7 +21,7 @@ def test_real_field_round_trip_preserves_shape_and_values() -> None:
restored = unvec(vector, SHAPE) restored = unvec(vector, SHAPE)
assert restored is not None assert restored is not None
assert restored.shape == (3, *SHAPE) assert restored.shape == (3, *SHAPE)
assert_array_equal(restored, FIELD) assert_close(restored, FIELD)
def test_complex_field_round_trip_preserves_shape_and_values() -> None: def test_complex_field_round_trip_preserves_shape_and_values() -> None:
@ -29,7 +30,7 @@ def test_complex_field_round_trip_preserves_shape_and_values() -> None:
restored = unvec(vector, SHAPE) restored = unvec(vector, SHAPE)
assert restored is not None assert restored is not None
assert restored.shape == (3, *SHAPE) assert restored.shape == (3, *SHAPE)
assert_allclose(restored, COMPLEX_FIELD) assert_close(restored, COMPLEX_FIELD)
def test_unvec_with_two_components_round_trips_vector() -> None: def test_unvec_with_two_components_round_trips_vector() -> None:
@ -37,9 +38,9 @@ def test_unvec_with_two_components_round_trips_vector() -> None:
field = unvec(vector, SHAPE, nvdim=2) field = unvec(vector, SHAPE, nvdim=2)
assert field is not None assert field is not None
assert field.shape == (2, *SHAPE) assert field.shape == (2, *SHAPE)
assert_array_equal(vec(field), vector) assert_close(vec(field), vector)
def test_vec_accepts_arraylike_input() -> None: def test_vec_accepts_arraylike_input() -> None:
arraylike = [[[1, 2], [3, 4]], [[5, 6], [7, 8]]] arraylike = [[[1, 2], [3, 4]], [[5, 6], [7, 8]]]
assert_array_equal(vec(arraylike), numpy.ravel(arraylike, order='C')) assert_close(vec(arraylike), numpy.ravel(arraylike, order='C'))

5
meanas/test/test_fdtd.py
View file

@ -7,7 +7,7 @@ from numpy.typing import NDArray
#from numpy.testing import assert_allclose, assert_array_equal #from numpy.testing import assert_allclose, assert_array_equal
from .. import fdtd from .. import fdtd
from .utils import assert_close, assert_fields_close, PRNG from .utils import assert_close, assert_fields_close, make_prng
from .conftest import FixtureRequest from .conftest import FixtureRequest
@ -179,13 +179,14 @@ def j_distribution(
shape: tuple[int, ...], shape: tuple[int, ...],
j_mag: float, j_mag: float,
) -> NDArray[numpy.float64]: ) -> NDArray[numpy.float64]:
prng = make_prng()
j = numpy.zeros(shape) j = numpy.zeros(shape)
if request.param == 'center': if request.param == 'center':
j[:, shape[1] // 2, shape[2] // 2, shape[3] // 2] = j_mag j[:, shape[1] // 2, shape[2] // 2, shape[3] // 2] = j_mag
elif request.param == '000': elif request.param == '000':
j[:, 0, 0, 0] = j_mag j[:, 0, 0, 0] = j_mag
elif request.param == 'random': elif request.param == 'random':
j[:] = PRNG.uniform(low=-j_mag, high=j_mag, size=shape) j[:] = prng.uniform(low=-j_mag, high=j_mag, size=shape)
return j return j

15
meanas/test/test_fdtd_base.py
View file

@ -1,14 +1,15 @@
import numpy import numpy
from numpy.testing import assert_allclose
from ..fdmath import functional as fd_functional from ..fdmath import functional as fd_functional
from ..fdtd import base from ..fdtd import base
from ._test_builders import real_ramp
from .utils import assert_close
DT = 0.25 DT = 0.25
SHAPE = (3, 2, 2, 2) SHAPE = (3, 2, 2, 2)
E_FIELD = numpy.arange(numpy.prod(SHAPE), dtype=float).reshape(SHAPE, order='C') / 5.0 E_FIELD = real_ramp(SHAPE, scale=1 / 5)
H_FIELD = (numpy.arange(numpy.prod(SHAPE), dtype=float).reshape(SHAPE, order='C') + 1.0) / 7.0 H_FIELD = real_ramp(SHAPE, scale=1 / 7, offset=1 / 7)
EPSILON = 1.5 + E_FIELD / 10.0 EPSILON = 1.5 + E_FIELD / 10.0
MU_FIELD = 2.0 + H_FIELD / 8.0 MU_FIELD = 2.0 + H_FIELD / 8.0
MU_SCALAR = 3.0 MU_SCALAR = 3.0
@ -20,7 +21,7 @@ def test_maxwell_e_without_dxes_matches_unit_spacing_update() -> None:
updated = updater(E_FIELD.copy(), H_FIELD.copy(), EPSILON) updated = updater(E_FIELD.copy(), H_FIELD.copy(), EPSILON)
assert_allclose(updated, expected) assert_close(updated, expected)
def test_maxwell_h_without_dxes_and_without_mu_matches_unit_spacing_update() -> None: def test_maxwell_h_without_dxes_and_without_mu_matches_unit_spacing_update() -> None:
@ -29,7 +30,7 @@ def test_maxwell_h_without_dxes_and_without_mu_matches_unit_spacing_update() ->
updated = updater(E_FIELD.copy(), H_FIELD.copy()) updated = updater(E_FIELD.copy(), H_FIELD.copy())
assert_allclose(updated, expected) assert_close(updated, expected)
def test_maxwell_h_without_dxes_accepts_scalar_and_field_mu() -> None: def test_maxwell_h_without_dxes_accepts_scalar_and_field_mu() -> None:
@ -37,8 +38,8 @@ def test_maxwell_h_without_dxes_accepts_scalar_and_field_mu() -> None:
updated_scalar = updater(E_FIELD.copy(), H_FIELD.copy(), MU_SCALAR) updated_scalar = updater(E_FIELD.copy(), H_FIELD.copy(), MU_SCALAR)
expected_scalar = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) / MU_SCALAR expected_scalar = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) / MU_SCALAR
assert_allclose(updated_scalar, expected_scalar) assert_close(updated_scalar, expected_scalar)
updated_field = updater(E_FIELD.copy(), H_FIELD.copy(), MU_FIELD) updated_field = updater(E_FIELD.copy(), H_FIELD.copy(), MU_FIELD)
expected_field = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) / MU_FIELD expected_field = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) / MU_FIELD
assert_allclose(updated_field, expected_field) assert_close(updated_field, expected_field)

62
meanas/test/test_fdtd_boundaries.py Normal file
View file

@ -0,0 +1,62 @@
import numpy
import pytest
from numpy.testing import assert_allclose
from ..fdtd.boundaries import conducting_boundary
def _axis_index(axis: int, index: int) -> tuple[slice | int, ...]:
coords: list[slice | int] = [slice(None), slice(None), slice(None)]
coords[axis] = index
return tuple(coords)
@pytest.mark.parametrize('direction', [0, 1, 2])
@pytest.mark.parametrize('polarity', [-1, 1])
def test_conducting_boundary_updates_expected_faces(direction: int, polarity: int) -> None:
e = numpy.arange(3 * 4 * 4 * 4, dtype=float).reshape(3, 4, 4, 4)
h = e.copy()
e0 = e.copy()
h0 = h.copy()
update_e, update_h = conducting_boundary(direction, polarity)
update_e(e)
update_h(h)
dirs = [0, 1, 2]
dirs.remove(direction)
u, v = dirs
if polarity < 0:
boundary = _axis_index(direction, 0)
shifted1 = _axis_index(direction, 1)
assert_allclose(e[direction][boundary], 0)
assert_allclose(e[u][boundary], e0[u][shifted1])
assert_allclose(e[v][boundary], e0[v][shifted1])
assert_allclose(h[direction][boundary], h0[direction][shifted1])
assert_allclose(h[u][boundary], 0)
assert_allclose(h[v][boundary], 0)
else:
boundary = _axis_index(direction, -1)
shifted1 = _axis_index(direction, -2)
shifted2 = _axis_index(direction, -3)
assert_allclose(e[direction][boundary], -e0[direction][shifted2])
assert_allclose(e[direction][shifted1], 0)
assert_allclose(e[u][boundary], e0[u][shifted1])
assert_allclose(e[v][boundary], e0[v][shifted1])
assert_allclose(h[direction][boundary], h0[direction][shifted1])
assert_allclose(h[u][boundary], -h0[u][shifted2])
assert_allclose(h[u][shifted1], 0)
assert_allclose(h[v][boundary], -h0[v][shifted2])
assert_allclose(h[v][shifted1], 0)
@pytest.mark.parametrize(
('direction', 'polarity'),
[(-1, 1), (3, 1), (0, 0)],
)
def test_conducting_boundary_rejects_invalid_arguments(direction: int, polarity: int) -> None:
with pytest.raises(Exception):
conducting_boundary(direction, polarity)

21
meanas/test/test_fdtd_energy.py
View file

@ -1,13 +1,14 @@
import numpy import numpy
from numpy.testing import assert_allclose
from .. import fdtd from .. import fdtd
from ..fdtd import energy as fdtd_energy from ..fdtd import energy as fdtd_energy
from ._test_builders import real_ramp, unit_dxes
from .utils import assert_close
SHAPE = (2, 2, 2) SHAPE = (2, 2, 2)
DT = 0.25 DT = 0.25
UNIT_DXES = tuple(tuple(numpy.ones(length) for length in SHAPE) for _ in range(2)) UNIT_DXES = unit_dxes(SHAPE)
DXES = ( DXES = (
( (
numpy.array([1.0, 1.5]), numpy.array([1.0, 1.5]),
@ -20,11 +21,11 @@ DXES = (
numpy.array([0.7, 1.3]), numpy.array([0.7, 1.3]),
), ),
) )
E0 = numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') E0 = real_ramp((3, *SHAPE))
E1 = E0 + 0.5 E1 = E0 + 0.5
E2 = E0 + 1.0 E2 = E0 + 1.0
E3 = E0 + 1.5 E3 = E0 + 1.5
H0 = (numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') + 2.0) / 3.0 H0 = real_ramp((3, *SHAPE), scale=1 / 3, offset=2 / 3)
H1 = H0 + 0.25 H1 = H0 + 0.25
H2 = H0 + 0.5 H2 = H0 + 0.5
H3 = H0 + 0.75 H3 = H0 + 0.75
@ -36,14 +37,14 @@ MU = 1.5 + H0 / 10.0
def test_poynting_default_spacing_matches_explicit_unit_spacing() -> None: def test_poynting_default_spacing_matches_explicit_unit_spacing() -> None:
default_spacing = fdtd.poynting(E1, H1) default_spacing = fdtd.poynting(E1, H1)
explicit_spacing = fdtd.poynting(E1, H1, dxes=UNIT_DXES) explicit_spacing = fdtd.poynting(E1, H1, dxes=UNIT_DXES)
assert_allclose(default_spacing, explicit_spacing) assert_close(default_spacing, explicit_spacing)
def test_poynting_divergence_matches_precomputed_poynting_vector() -> None: def test_poynting_divergence_matches_precomputed_poynting_vector() -> None:
s = fdtd.poynting(E2, H2, dxes=DXES) s = fdtd.poynting(E2, H2, dxes=DXES)
from_fields = fdtd.poynting_divergence(e=E2, h=H2, dxes=DXES) from_fields = fdtd.poynting_divergence(e=E2, h=H2, dxes=DXES)
from_vector = fdtd.poynting_divergence(s=s) from_vector = fdtd.poynting_divergence(s=s)
assert_allclose(from_fields, from_vector) assert_close(from_fields, from_vector)
def test_delta_energy_h2e_matches_direct_dxmul_formula() -> None: def test_delta_energy_h2e_matches_direct_dxmul_formula() -> None:
@ -64,7 +65,7 @@ def test_delta_energy_h2e_matches_direct_dxmul_formula() -> None:
mu=MU, mu=MU,
dxes=DXES, dxes=DXES,
) )
assert_allclose(actual, expected) assert_close(actual, expected)
def test_delta_energy_e2h_matches_direct_dxmul_formula() -> None: def test_delta_energy_e2h_matches_direct_dxmul_formula() -> None:
@ -85,14 +86,14 @@ def test_delta_energy_e2h_matches_direct_dxmul_formula() -> None:
mu=MU, mu=MU,
dxes=DXES, dxes=DXES,
) )
assert_allclose(actual, expected) assert_close(actual, expected)
def test_delta_energy_j_defaults_to_unit_cell_volume() -> None: def test_delta_energy_j_defaults_to_unit_cell_volume() -> None:
expected = (J0 * E1).sum(axis=0) expected = (J0 * E1).sum(axis=0)
assert_allclose(fdtd.delta_energy_j(j0=J0, e1=E1), expected) assert_close(fdtd.delta_energy_j(j0=J0, e1=E1), expected)
def test_dxmul_defaults_to_unit_materials_and_spacing() -> None: def test_dxmul_defaults_to_unit_materials_and_spacing() -> None:
expected = E1.sum(axis=0) + H1.sum(axis=0) expected = E1.sum(axis=0) + H1.sum(axis=0)
assert_allclose(fdtd_energy.dxmul(E1, H1), expected) assert_close(fdtd_energy.dxmul(E1, H1), expected)

42
meanas/test/test_fdtd_misc.py Normal file
View file

@ -0,0 +1,42 @@
import numpy
import pytest
from ..fdtd.misc import gaussian_beam, gaussian_packet, ricker_pulse
@pytest.mark.parametrize('one_sided', [False, True])
def test_gaussian_packet_accepts_array_input(one_sided: bool) -> None:
dt = 0.01
source, delay = gaussian_packet(1.55, 0.1, dt, one_sided=one_sided)
steps = numpy.array([0, int(numpy.ceil(delay / dt)) + 5])
envelope, cc, ss = source(steps)
assert envelope.shape == (2,)
assert numpy.isfinite(envelope).all()
assert numpy.isfinite(cc).all()
assert numpy.isfinite(ss).all()
if one_sided:
assert envelope[-1] == pytest.approx(1.0)
def test_ricker_pulse_returns_finite_values() -> None:
source, delay = ricker_pulse(1.55, 0.01)
envelope, cc, ss = source(numpy.array([0, 1, 2]))
assert numpy.isfinite(delay)
assert numpy.isfinite(envelope).all()
assert numpy.isfinite(cc).all()
assert numpy.isfinite(ss).all()
def test_gaussian_beam_centered_grid_is_finite_and_normalized() -> None:
beam = gaussian_beam(
xyz=[numpy.linspace(-1, 1, 3), numpy.linspace(-1, 1, 3), numpy.linspace(-1, 1, 3)],
center=[0, 0, 0],
waist_radius=1.0,
wl=1.55,
)
row = beam[:, :, beam.shape[2] // 2]
assert numpy.isfinite(beam).all()
assert numpy.linalg.norm(row) == pytest.approx(1.0)

44
meanas/test/test_fdtd_pml.py Normal file
View file

@ -0,0 +1,44 @@
import numpy
import pytest
from ..fdtd.pml import cpml_params, updates_with_cpml
@pytest.mark.parametrize(
('axis', 'polarity', 'thickness', 'epsilon_eff'),
[(3, 1, 4, 1.0), (0, 0, 4, 1.0), (0, 1, 2, 1.0), (0, 1, 4, 0.0)],
)
def test_cpml_params_reject_invalid_arguments(axis: int, polarity: int, thickness: int, epsilon_eff: float) -> None:
with pytest.raises(Exception):
cpml_params(axis=axis, polarity=polarity, dt=0.1, thickness=thickness, epsilon_eff=epsilon_eff)
def test_cpml_params_shapes_and_region() -> None:
params = cpml_params(axis=1, polarity=1, dt=0.1, thickness=3)
p0e, p1e, p2e = params['param_e']
p0h, p1h, p2h = params['param_h']
assert p0e.shape == (1, 3, 1)
assert p1e.shape == (1, 3, 1)
assert p2e.shape == (1, 3, 1)
assert p0h.shape == (1, 3, 1)
assert p1h.shape == (1, 3, 1)
assert p2h.shape == (1, 3, 1)
assert params['region'][1] == slice(-3, None)
def test_updates_with_cpml_keeps_zero_fields_zero() -> None:
shape = (3, 4, 4, 4)
epsilon = numpy.ones(shape, dtype=float)
e = numpy.zeros(shape, dtype=float)
h = numpy.zeros(shape, dtype=float)
dxes = [[numpy.ones(4), numpy.ones(4), numpy.ones(4)] for _ in range(2)]
params = [[None, None] for _ in range(3)]
params[0][0] = cpml_params(axis=0, polarity=-1, dt=0.1, thickness=3)
update_e, update_h = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
update_e(e, h, epsilon)
update_h(e, h)
assert not e.any()
assert not h.any()

20
meanas/test/test_import_fallbacks.py
View file

@ -4,6 +4,16 @@ import pathlib
import meanas import meanas
from ..fdfd import bloch from ..fdfd import bloch
from .utils import assert_close
def _reload(module):
return importlib.reload(module)
def _restore_reloaded(monkeypatch, module):
monkeypatch.undo()
return _reload(module)
def test_meanas_import_survives_readme_open_failure(monkeypatch) -> None: # type: ignore[no-untyped-def] def test_meanas_import_survives_readme_open_failure(monkeypatch) -> None: # type: ignore[no-untyped-def]
@ -15,14 +25,13 @@ def test_meanas_import_survives_readme_open_failure(monkeypatch) -> None: # typ
return original_open(self, *args, **kwargs) return original_open(self, *args, **kwargs)
monkeypatch.setattr(pathlib.Path, 'open', failing_open) monkeypatch.setattr(pathlib.Path, 'open', failing_open)
reloaded = importlib.reload(meanas) reloaded = _reload(meanas)
assert reloaded.__version__ == '0.10' assert reloaded.__version__ == '0.10'
assert reloaded.__author__ == 'Jan Petykiewicz' assert reloaded.__author__ == 'Jan Petykiewicz'
assert reloaded.__doc__ is not None assert reloaded.__doc__ is not None
monkeypatch.undo() _restore_reloaded(monkeypatch, meanas)
importlib.reload(meanas)
def test_bloch_reloads_with_numpy_fft_when_pyfftw_is_unavailable(monkeypatch) -> None: # type: ignore[no-untyped-def] def test_bloch_reloads_with_numpy_fft_when_pyfftw_is_unavailable(monkeypatch) -> None: # type: ignore[no-untyped-def]
@ -34,10 +43,9 @@ def test_bloch_reloads_with_numpy_fft_when_pyfftw_is_unavailable(monkeypatch) ->
return original_import(name, globals, locals, fromlist, level) return original_import(name, globals, locals, fromlist, level)
monkeypatch.setattr(builtins, '__import__', fake_import) monkeypatch.setattr(builtins, '__import__', fake_import)
reloaded = importlib.reload(bloch) reloaded = _reload(bloch)
assert reloaded.fftn.__module__ == 'numpy.fft' assert reloaded.fftn.__module__ == 'numpy.fft'
assert reloaded.ifftn.__module__ == 'numpy.fft' assert reloaded.ifftn.__module__ == 'numpy.fft'
monkeypatch.undo() _restore_reloaded(monkeypatch, bloch)
importlib.reload(bloch)

238
meanas/test/test_regressions.py
View file

@ -1,238 +0,0 @@
import numpy
import pytest # type: ignore
from numpy.testing import assert_allclose
from scipy import sparse
from ..eigensolvers import power_iteration, rayleigh_quotient_iteration, signed_eigensolve
from ..fdfd import eme, farfield
from ..fdtd.boundaries import conducting_boundary
from ..fdtd.misc import gaussian_beam, gaussian_packet, ricker_pulse
from ..fdtd.pml import cpml_params, updates_with_cpml
def _axis_index(axis: int, index: int) -> tuple[slice | int, ...]:
coords: list[slice | int] = [slice(None), slice(None), slice(None)]
coords[axis] = index
return tuple(coords)
@pytest.mark.parametrize('one_sided', [False, True])
def test_gaussian_packet_accepts_array_input(one_sided: bool) -> None:
dt = 0.01
source, delay = gaussian_packet(1.55, 0.1, dt, one_sided=one_sided)
steps = numpy.array([0, int(numpy.ceil(delay / dt)) + 5])
envelope, cc, ss = source(steps)
assert envelope.shape == (2,)
assert numpy.isfinite(envelope).all()
assert numpy.isfinite(cc).all()
assert numpy.isfinite(ss).all()
if one_sided:
assert envelope[-1] == pytest.approx(1.0)
def test_ricker_pulse_returns_finite_values() -> None:
source, delay = ricker_pulse(1.55, 0.01)
envelope, cc, ss = source(numpy.array([0, 1, 2]))
assert numpy.isfinite(delay)
assert numpy.isfinite(envelope).all()
assert numpy.isfinite(cc).all()
assert numpy.isfinite(ss).all()
def test_gaussian_beam_centered_grid_is_finite_and_normalized() -> None:
beam = gaussian_beam(
xyz=[numpy.linspace(-1, 1, 3), numpy.linspace(-1, 1, 3), numpy.linspace(-1, 1, 3)],
center=[0, 0, 0],
waist_radius=1.0,
wl=1.55,
)
row = beam[:, :, beam.shape[2] // 2]
assert numpy.isfinite(beam).all()
assert numpy.linalg.norm(row) == pytest.approx(1.0)
@pytest.mark.parametrize('direction', [0, 1, 2])
@pytest.mark.parametrize('polarity', [-1, 1])
def test_conducting_boundary_updates_expected_faces(direction: int, polarity: int) -> None:
e = numpy.arange(3 * 4 * 4 * 4, dtype=float).reshape(3, 4, 4, 4)
h = e.copy()
e0 = e.copy()
h0 = h.copy()
update_e, update_h = conducting_boundary(direction, polarity)
update_e(e)
update_h(h)
dirs = [0, 1, 2]
dirs.remove(direction)
u, v = dirs
if polarity < 0:
boundary = _axis_index(direction, 0)
shifted1 = _axis_index(direction, 1)
assert_allclose(e[direction][boundary], 0)
assert_allclose(e[u][boundary], e0[u][shifted1])
assert_allclose(e[v][boundary], e0[v][shifted1])
assert_allclose(h[direction][boundary], h0[direction][shifted1])
assert_allclose(h[u][boundary], 0)
assert_allclose(h[v][boundary], 0)
else:
boundary = _axis_index(direction, -1)
shifted1 = _axis_index(direction, -2)
shifted2 = _axis_index(direction, -3)
assert_allclose(e[direction][boundary], -e0[direction][shifted2])
assert_allclose(e[direction][shifted1], 0)
assert_allclose(e[u][boundary], e0[u][shifted1])
assert_allclose(e[v][boundary], e0[v][shifted1])
assert_allclose(h[direction][boundary], h0[direction][shifted1])
assert_allclose(h[u][boundary], -h0[u][shifted2])
assert_allclose(h[u][shifted1], 0)
assert_allclose(h[v][boundary], -h0[v][shifted2])
assert_allclose(h[v][shifted1], 0)
@pytest.mark.parametrize(
('direction', 'polarity'),
[(-1, 1), (3, 1), (0, 0)],
)
def test_conducting_boundary_rejects_invalid_arguments(direction: int, polarity: int) -> None:
with pytest.raises(Exception):
conducting_boundary(direction, polarity)
def test_get_abcd_returns_sparse_block_matrix(monkeypatch: pytest.MonkeyPatch) -> None:
t = numpy.array([[1.0, 0.0], [0.0, 2.0]])
r = numpy.array([[0.0, 0.1], [0.2, 0.0]])
monkeypatch.setattr(eme, 'get_tr', lambda *args, **kwargs: (t, r))
abcd = eme.get_abcd(None, None, None, None)
assert sparse.issparse(abcd)
assert abcd.shape == (4, 4)
assert numpy.isfinite(abcd.toarray()).all()
def test_get_s_force_reciprocal_symmetrizes_output(monkeypatch: pytest.MonkeyPatch) -> None:
left = object()
right = object()
def fake_get_tr(_eh_l, wavenumbers_l, _eh_r, _wavenumbers_r, **kwargs):
if wavenumbers_l is left:
return (
numpy.array([[1.0, 2.0], [0.5, 1.0]]),
numpy.array([[0.0, 1.0], [2.0, 0.0]]),
)
return (
numpy.array([[1.0, -1.0], [0.0, 1.0]]),
numpy.array([[0.0, 0.5], [1.5, 0.0]]),
)
monkeypatch.setattr(eme, 'get_tr', fake_get_tr)
ss = eme.get_s(None, left, None, right, force_reciprocal=True)
assert_allclose(ss, ss.T)
def test_farfield_roundtrip_supports_rectangular_arrays() -> None:
e_near = [numpy.zeros((4, 8), dtype=complex), numpy.zeros((4, 8), dtype=complex)]
h_near = [numpy.zeros((4, 8), dtype=complex), numpy.zeros((4, 8), dtype=complex)]
e_near[0][1, 3] = 1.0 + 0.25j
h_near[1][2, 5] = -0.5j
ff = farfield.near_to_farfield(e_near, h_near, dx=0.2, dy=0.3, padded_size=(4, 8))
restored = farfield.far_to_nearfield(
[field.copy() for field in ff['E']],
[field.copy() for field in ff['H']],
ff['dkx'],
ff['dky'],
padded_size=(4, 8),
)
assert isinstance(ff['dkx'], float)
assert isinstance(ff['dky'], float)
assert ff['E'][0].shape == (4, 8)
assert restored['E'][0].shape == (4, 8)
assert restored['H'][0].shape == (4, 8)
assert restored['dx'] == pytest.approx(0.2)
assert restored['dy'] == pytest.approx(0.3)
assert numpy.isfinite(restored['E'][0]).all()
assert numpy.isfinite(restored['H'][0]).all()
@pytest.mark.parametrize(
('axis', 'polarity', 'thickness', 'epsilon_eff'),
[(3, 1, 4, 1.0), (0, 0, 4, 1.0), (0, 1, 2, 1.0), (0, 1, 4, 0.0)],
)
def test_cpml_params_reject_invalid_arguments(axis: int, polarity: int, thickness: int, epsilon_eff: float) -> None:
with pytest.raises(Exception):
cpml_params(axis=axis, polarity=polarity, dt=0.1, thickness=thickness, epsilon_eff=epsilon_eff)
def test_cpml_params_shapes_and_region() -> None:
params = cpml_params(axis=1, polarity=1, dt=0.1, thickness=3)
p0e, p1e, p2e = params['param_e']
p0h, p1h, p2h = params['param_h']
assert p0e.shape == (1, 3, 1)
assert p1e.shape == (1, 3, 1)
assert p2e.shape == (1, 3, 1)
assert p0h.shape == (1, 3, 1)
assert p1h.shape == (1, 3, 1)
assert p2h.shape == (1, 3, 1)
assert params['region'][1] == slice(-3, None)
def test_updates_with_cpml_keeps_zero_fields_zero() -> None:
shape = (3, 4, 4, 4)
epsilon = numpy.ones(shape, dtype=float)
e = numpy.zeros(shape, dtype=float)
h = numpy.zeros(shape, dtype=float)
dxes = [[numpy.ones(4), numpy.ones(4), numpy.ones(4)] for _ in range(2)]
params = [[None, None] for _ in range(3)]
params[0][0] = cpml_params(axis=0, polarity=-1, dt=0.1, thickness=3)
update_e, update_h = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
update_e(e, h, epsilon)
update_h(e, h)
assert not e.any()
assert not h.any()
def test_power_iteration_finds_dominant_mode() -> None:
operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
eigval, eigvec = power_iteration(operator, guess_vector=numpy.ones(4, dtype=complex), iterations=20)
assert eigval == pytest.approx(5.0, rel=1e-6)
assert abs(eigvec[0]) > abs(eigvec[1])
def test_rayleigh_quotient_iteration_refines_known_mode() -> None:
operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
def solver(matrix: sparse.spmatrix, rhs: numpy.ndarray) -> numpy.ndarray:
return numpy.linalg.lstsq(matrix.toarray(), rhs, rcond=None)[0]
eigval, eigvec = rayleigh_quotient_iteration(
operator,
numpy.array([1.0, 0.1, 0.0, 0.0], dtype=complex),
iterations=8,
solver=solver,
)
residual = numpy.linalg.norm(operator @ eigvec - eigval * eigvec)
assert eigval == pytest.approx(3.0, rel=1e-6)
assert residual < 1e-8
def test_signed_eigensolve_returns_largest_positive_modes() -> None:
operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
eigvals, eigvecs = signed_eigensolve(operator, how_many=2)
assert_allclose(eigvals, [3.0, 5.0], atol=1e-6)
assert eigvecs.shape == (4, 2)

4
meanas/test/utils.py
View file

@ -2,7 +2,8 @@ import numpy
from numpy.typing import NDArray from numpy.typing import NDArray
PRNG = numpy.random.RandomState(12345) def make_prng(seed: int = 12345) -> numpy.random.RandomState:
return numpy.random.RandomState(seed)
def assert_fields_close( def assert_fields_close(
@ -29,4 +30,3 @@ def assert_close(
**kwargs, **kwargs,
) -> None: ) -> None:
numpy.testing.assert_allclose(x, y, *args, **kwargs) numpy.testing.assert_allclose(x, y, *args, **kwargs)