From 593098bf8fc92b9ee8920a2783c08393e516bddf Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Fri, 17 Apr 2026 20:30:28 -0700 Subject: [PATCH] [fdfd.functional] fix handling of mu in e_full and m2j sign --- meanas/fdfd/functional.py | 12 +-- meanas/test/test_fdfd_functional.py | 130 ++++++++++++++++++++++++++++ 2 files changed, 136 insertions(+), 6 deletions(-) create mode 100644 meanas/test/test_fdfd_functional.py diff --git a/meanas/fdfd/functional.py b/meanas/fdfd/functional.py index 440daf2..9d98798 100644 --- a/meanas/fdfd/functional.py +++ b/meanas/fdfd/functional.py @@ -41,8 +41,8 @@ def e_full( curls = ch(ce(e)) return cfdfield_t(curls - omega ** 2 * epsilon * e) - def op_mu(e: cfdfield) -> cfdfield_t: - curls = ch(mu * ce(e)) # type: ignore # mu = None ok because we don't return the function + def op_mu(e: cfdfield_t) -> cfdfield_t: + curls = ch(ce(e) / mu) # type: ignore # mu = None ok because we don't return the function return cfdfield_t(curls - omega ** 2 * epsilon * e) if mu is None: @@ -138,12 +138,12 @@ def m2j( """ ch = curl_back(dxes[1]) - def m2j_mu(m: cfdfield) -> cfdfield_t: - J = ch(m / mu) / (-1j * omega) # type: ignore # mu=None ok + def m2j_mu(m: cfdfield_t) -> cfdfield_t: + J = ch(m / mu) / (1j * omega) # type: ignore # mu=None ok return cfdfield_t(J) - def m2j_1(m: cfdfield) -> cfdfield_t: - J = ch(m) / (-1j * omega) + def m2j_1(m: cfdfield_t) -> cfdfield_t: + J = ch(m) / (1j * omega) return cfdfield_t(J) if mu is None: diff --git a/meanas/test/test_fdfd_functional.py b/meanas/test/test_fdfd_functional.py new file mode 100644 index 0000000..f4fd4bb --- /dev/null +++ b/meanas/test/test_fdfd_functional.py @@ -0,0 +1,130 @@ +import numpy +from numpy.testing import assert_allclose + +from ..fdmath import vec, unvec +from ..fdfd import functional, operators + + +OMEGA = 1 / 1500 +SHAPE = (2, 3, 2) +ATOL = 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: + assert_allclose(actual, expected, atol=ATOL, rtol=RTOL) + + +def test_e_full_matches_sparse_operator_without_mu() -> None: + matrix_result = apply_matrix( + operators.e_full(OMEGA, DXES, vec(EPSILON)), + E_FIELD, + ) + functional_result = functional.e_full(OMEGA, DXES, EPSILON)(E_FIELD) + + assert_fields_match(functional_result, matrix_result) + + +def test_e_full_matches_sparse_operator_with_mu() -> None: + matrix_result = apply_matrix( + operators.e_full(OMEGA, DXES, vec(EPSILON), vec(MU)), + E_FIELD, + ) + functional_result = functional.e_full(OMEGA, DXES, EPSILON, MU)(E_FIELD) + + assert_fields_match(functional_result, matrix_result) + + +def test_eh_full_matches_sparse_operator_with_mu() -> None: + matrix_result = operators.eh_full(OMEGA, DXES, vec(EPSILON), vec(MU)) @ 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, MU)(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)), + E_FIELD, + ) + functional_result = functional.e2h(OMEGA, DXES, MU)(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), + H_FIELD, + ) + functional_result = functional.m2j(OMEGA, DXES)(H_FIELD) + + assert_fields_match(functional_result, matrix_result) + + +def test_m2j_matches_sparse_operator_with_mu() -> None: + matrix_result = apply_matrix( + operators.m2j(OMEGA, DXES, vec(MU)), + H_FIELD, + ) + functional_result = functional.m2j(OMEGA, DXES, MU)(H_FIELD) + + assert_fields_match(functional_result, matrix_result) + + +def test_e_tfsf_source_matches_sparse_operator_without_mu() -> None: + matrix_result = apply_matrix( + operators.e_tfsf_source(vec(TF_REGION), OMEGA, DXES, vec(EPSILON)), + E_FIELD, + ) + functional_result = functional.e_tfsf_source(TF_REGION, OMEGA, DXES, EPSILON)(E_FIELD) + + assert_fields_match(functional_result, matrix_result) + + +def test_e_tfsf_source_matches_sparse_operator_with_mu() -> None: + matrix_result = apply_matrix( + operators.e_tfsf_source(vec(TF_REGION), OMEGA, DXES, vec(EPSILON), vec(MU)), + E_FIELD, + ) + functional_result = functional.e_tfsf_source(TF_REGION, OMEGA, DXES, EPSILON, MU)(E_FIELD) + + assert_fields_match(functional_result, matrix_result) + + +def test_poynting_e_cross_h_matches_sparse_operator() -> None: + matrix_result = apply_matrix( + operators.poynting_e_cross(vec(E_FIELD), DXES), + H_FIELD, + ) + functional_result = functional.poynting_e_cross_h(DXES)(E_FIELD, H_FIELD) + + assert_fields_match(functional_result, matrix_result)