[fdfd.operators] fix eh_full for non-None mu
This commit is contained in:
parent
0ff23542ac
commit
0afe2297b0
4 changed files with 178 additions and 4 deletions
|
|
@ -236,10 +236,12 @@ def eh_full(
|
|||
else:
|
||||
pm = sparse.diags_array(numpy.where(pmc, 0, 1)) # set pm to (not PMC)
|
||||
|
||||
iwe = pe @ (1j * omega * sparse.diags_array(epsilon)) @ pe
|
||||
iwm = 1j * omega
|
||||
if mu is not None:
|
||||
iwm *= sparse.diags_array(mu)
|
||||
iwe = pe @ (1j * omega * sparse.diags(epsilon)) @ pe
|
||||
if mu is None:
|
||||
iwm = 1j * omega * sparse.eye(epsilon.size)
|
||||
else:
|
||||
iwm = 1j * omega * sparse.diags(mu)
|
||||
|
||||
iwm = pm @ iwm @ pm
|
||||
|
||||
A1 = pe @ curl_back(dxes[1]) @ pm
|
||||
|
|
|
|||
|
|
@ -49,6 +49,21 @@ def _apply_matrix(op: operators.sparse.spmatrix, field: numpy.ndarray, shape: tu
|
|||
return unvec(op @ vec(field), shape)
|
||||
|
||||
|
||||
def _dense_e_full(mu: numpy.ndarray | None) -> numpy.ndarray:
|
||||
ce = fd_functional.curl_forward(DXES[0])
|
||||
ch = fd_functional.curl_back(DXES[1])
|
||||
pe = numpy.where(PEC, 0.0, 1.0)
|
||||
pm = numpy.where(PMC, 0.0, 1.0)
|
||||
|
||||
masked_e = pe * E_FIELD
|
||||
curl_term = ce(masked_e)
|
||||
if mu is not None:
|
||||
curl_term = curl_term / mu
|
||||
curl_term = pm * curl_term
|
||||
curl_term = ch(curl_term)
|
||||
return pe * (curl_term - OMEGA**2 * EPSILON * masked_e)
|
||||
|
||||
|
||||
def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray:
|
||||
ce = fd_functional.curl_forward(DXES[0])
|
||||
ch = fd_functional.curl_back(DXES[1])
|
||||
|
|
@ -78,6 +93,55 @@ def test_h_full_matches_dense_reference_with_and_without_mu() -> None:
|
|||
assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
|
||||
|
||||
|
||||
def test_e_full_matches_dense_reference_with_masks() -> None:
|
||||
for mu in (None, MU):
|
||||
matrix_result = _apply_matrix(
|
||||
operators.e_full(OMEGA, DXES, vec(EPSILON), None if mu is None else vec(mu), vec(PEC), vec(PMC)),
|
||||
E_FIELD,
|
||||
SHAPE,
|
||||
)
|
||||
dense_result = _dense_e_full(mu)
|
||||
assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
|
||||
|
||||
|
||||
def test_h_full_without_masks_matches_dense_reference() -> None:
|
||||
ce = fd_functional.curl_forward(DXES[0])
|
||||
ch = fd_functional.curl_back(DXES[1])
|
||||
dense_result = ce(ch(H_FIELD) / EPSILON) - OMEGA**2 * MU * H_FIELD
|
||||
matrix_result = _apply_matrix(
|
||||
operators.h_full(OMEGA, DXES, vec(EPSILON), vec(MU)),
|
||||
H_FIELD,
|
||||
SHAPE,
|
||||
)
|
||||
assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
|
||||
|
||||
|
||||
def test_eh_full_matches_manual_block_operator_with_masks() -> None:
|
||||
pe = numpy.where(PEC, 0.0, 1.0)
|
||||
pm = numpy.where(PMC, 0.0, 1.0)
|
||||
ce = fd_functional.curl_forward(DXES[0])
|
||||
ch = fd_functional.curl_back(DXES[1])
|
||||
|
||||
matrix_result = operators.eh_full(OMEGA, DXES, vec(EPSILON), vec(MU), vec(PEC), vec(PMC)) @ numpy.concatenate(
|
||||
[vec(E_FIELD), vec(H_FIELD)],
|
||||
)
|
||||
matrix_e, matrix_h = (unvec(part, SHAPE) for part in numpy.split(matrix_result, 2))
|
||||
|
||||
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))
|
||||
|
||||
assert_allclose(matrix_e, dense_e, atol=1e-10, rtol=1e-10)
|
||||
assert_allclose(matrix_h, dense_h, atol=1e-10, rtol=1e-10)
|
||||
|
||||
|
||||
def test_e2h_pmc_mask_matches_masked_unmasked_result() -> None:
|
||||
pmc_complement = numpy.where(PMC, 0.0, 1.0)
|
||||
unmasked = _apply_matrix(operators.e2h(OMEGA, DXES, vec(MU)), E_FIELD, SHAPE)
|
||||
masked = _apply_matrix(operators.e2h(OMEGA, DXES, vec(MU), vec(PMC)), E_FIELD, SHAPE)
|
||||
|
||||
assert_allclose(masked, pmc_complement * unmasked, atol=1e-10, rtol=1e-10)
|
||||
|
||||
|
||||
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)
|
||||
e_cross_h = _apply_matrix(operators.poynting_e_cross(vec(E_FIELD), DXES), H_FIELD, SHAPE)
|
||||
|
|
|
|||
|
|
@ -70,6 +70,15 @@ def test_eh_full_matches_sparse_operator_with_mu() -> None:
|
|||
assert_fields_match(functional_h, matrix_h)
|
||||
|
||||
|
||||
def test_eh_full_matches_sparse_operator_without_mu() -> None:
|
||||
matrix_result = operators.eh_full(OMEGA, DXES, vec(EPSILON)) @ numpy.concatenate([vec(E_FIELD), vec(H_FIELD)])
|
||||
matrix_e, matrix_h = (unvec(part, SHAPE) for part in numpy.split(matrix_result, 2))
|
||||
functional_e, functional_h = functional.eh_full(OMEGA, DXES, EPSILON)(E_FIELD, H_FIELD)
|
||||
|
||||
assert_fields_match(functional_e, matrix_e)
|
||||
assert_fields_match(functional_h, matrix_h)
|
||||
|
||||
|
||||
def test_e2h_matches_sparse_operator_with_mu() -> None:
|
||||
matrix_result = apply_matrix(
|
||||
operators.e2h(OMEGA, DXES, vec(MU)),
|
||||
|
|
@ -80,6 +89,16 @@ def test_e2h_matches_sparse_operator_with_mu() -> None:
|
|||
assert_fields_match(functional_result, matrix_result)
|
||||
|
||||
|
||||
def test_e2h_matches_sparse_operator_without_mu() -> None:
|
||||
matrix_result = apply_matrix(
|
||||
operators.e2h(OMEGA, DXES),
|
||||
E_FIELD,
|
||||
)
|
||||
functional_result = functional.e2h(OMEGA, DXES)(E_FIELD)
|
||||
|
||||
assert_fields_match(functional_result, matrix_result)
|
||||
|
||||
|
||||
def test_m2j_matches_sparse_operator_without_mu() -> None:
|
||||
matrix_result = apply_matrix(
|
||||
operators.m2j(OMEGA, DXES),
|
||||
|
|
|
|||
89
meanas/test/test_fdmath_operators.py
Normal file
89
meanas/test/test_fdmath_operators.py
Normal file
|
|
@ -0,0 +1,89 @@
|
|||
import numpy
|
||||
import pytest
|
||||
from numpy.testing import assert_allclose, assert_array_equal
|
||||
|
||||
from ..fdmath import operators, unvec, vec
|
||||
|
||||
|
||||
SHAPE = (2, 3, 2)
|
||||
SCALAR_FIELD = numpy.arange(numpy.prod(SHAPE), dtype=float).reshape(SHAPE, order='C')
|
||||
VECTOR_LEFT = (numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') + 0.5)
|
||||
VECTOR_RIGHT = (2.0 + numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') / 3.0)
|
||||
|
||||
|
||||
def _apply_scalar_matrix(op: operators.sparse.spmatrix) -> numpy.ndarray:
|
||||
return (op @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C')
|
||||
|
||||
|
||||
def _mirrored_indices(size: int, shift_distance: int) -> numpy.ndarray:
|
||||
indices = numpy.arange(size) + shift_distance
|
||||
indices = numpy.where(indices >= size, 2 * size - indices - 1, indices)
|
||||
indices = numpy.where(indices < 0, -1 - indices, indices)
|
||||
return indices
|
||||
|
||||
|
||||
@pytest.mark.parametrize(('axis', 'shift_distance'), [(0, 1), (1, -1), (2, 1)])
|
||||
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))
|
||||
expected = numpy.roll(SCALAR_FIELD, -shift_distance, axis=axis)
|
||||
assert_array_equal(matrix_result, expected)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(('axis', 'shift_distance'), [(0, 1), (1, -1), (2, 1)])
|
||||
def test_shift_with_mirror_matches_explicit_mirrored_indices(axis: int, shift_distance: int) -> None:
|
||||
matrix_result = _apply_scalar_matrix(operators.shift_with_mirror(axis, SHAPE, shift_distance))
|
||||
indices = [numpy.arange(length) for length in SHAPE]
|
||||
indices[axis] = _mirrored_indices(SHAPE[axis], shift_distance)
|
||||
expected = SCALAR_FIELD[numpy.ix_(*indices)]
|
||||
assert_array_equal(matrix_result, expected)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
('args', 'message'),
|
||||
[
|
||||
((0, (2,), 1), 'Invalid shape'),
|
||||
((3, SHAPE, 1), 'Invalid direction'),
|
||||
],
|
||||
)
|
||||
def test_shift_circ_rejects_invalid_arguments(args: tuple[int, tuple[int, ...], int], message: str) -> None:
|
||||
with pytest.raises(Exception, match=message):
|
||||
operators.shift_circ(*args)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
('args', 'message'),
|
||||
[
|
||||
((0, (2,), 1), 'Invalid shape'),
|
||||
((3, SHAPE, 1), 'Invalid direction'),
|
||||
((0, SHAPE, SHAPE[0]), 'too large'),
|
||||
],
|
||||
)
|
||||
def test_shift_with_mirror_rejects_invalid_arguments(args: tuple[int, tuple[int, ...], int], message: str) -> None:
|
||||
with pytest.raises(Exception, match=message):
|
||||
operators.shift_with_mirror(*args)
|
||||
|
||||
|
||||
def test_vec_cross_matches_pointwise_cross_product() -> None:
|
||||
matrix_result = unvec(operators.vec_cross(vec(VECTOR_LEFT)) @ vec(VECTOR_RIGHT), SHAPE)
|
||||
expected = numpy.empty_like(VECTOR_LEFT)
|
||||
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[2] = VECTOR_LEFT[0] * VECTOR_RIGHT[1] - VECTOR_LEFT[1] * VECTOR_RIGHT[0]
|
||||
assert_allclose(matrix_result, expected)
|
||||
|
||||
|
||||
def test_avg_forward_matches_half_sum_with_forward_neighbor() -> None:
|
||||
matrix_result = _apply_scalar_matrix(operators.avg_forward(1, SHAPE))
|
||||
expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, -1, axis=1))
|
||||
assert_allclose(matrix_result, expected)
|
||||
|
||||
|
||||
def test_avg_back_matches_half_sum_with_backward_neighbor() -> None:
|
||||
matrix_result = _apply_scalar_matrix(operators.avg_back(1, SHAPE))
|
||||
expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, 1, axis=1))
|
||||
assert_allclose(matrix_result, expected)
|
||||
|
||||
|
||||
def test_avg_forward_rejects_invalid_shape() -> None:
|
||||
with pytest.raises(Exception, match='Invalid shape'):
|
||||
operators.avg_forward(0, (2,))
|
||||
Loading…
Add table
Add a link
Reference in a new issue