[fdfd.functional] fix handling of mu in e_full and m2j sign

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:
+    def op_mu(e: cfdfield_t) -> cfdfield_t:
-        curls = ch(mu * ce(e))          # type: ignore   # mu = None ok because we don't return the function
+        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:
+    def m2j_mu(m: cfdfield_t) -> cfdfield_t:
-        J = ch(m / mu) / (-1j * omega)          # type: ignore  # mu=None ok
+        J = ch(m / mu) / (1j * omega)           # type: ignore  # mu=None ok
         return cfdfield_t(J)
 
-    def m2j_1(m: cfdfield) -> cfdfield_t:
+    def m2j_1(m: cfdfield_t) -> cfdfield_t:
-        J = ch(m) / (-1j * omega)
+        J = ch(m) / (1j * omega)
         return cfdfield_t(J)
 
     if mu is None:

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)