Merge pull request #28 from JakobEliasWagner/feature/dolfinx-function…

…-interpolation

Feature dolfinx function interpolation

src/nos/data/helmholtz/solver/__init__.py

@@ -1,5 +1,7 @@
+from .adiabatic_absorber import AdiabaticAbsorber
 from .util import get_mesh
 
 __all__ = [
     "get_mesh",
+    "AdiabaticAbsorber",
 ]

src/nos/data/helmholtz/solver/adiabatic_absorber.py

@@ -0,0 +1,39 @@
+import numpy as np
+
+from nos.data.helmholtz.domain_properties import Description
+from nos.utility import BoundingBox2D
+
+from .wave_number_function import WaveNumberFunction
+
+
+class AdiabaticAbsorber(WaveNumberFunction):
+    """Adiabatic layer class to truncate a domain with trivial reflections.
+
+    Shown in Oskooi et al. "The failure of perfectly matched layers, and towards their redemption by adiabatic
+    absorbers" Adiabatic layers asymptotically approach an adiabatic limit of zero reflections. This is done by
+    gradually increasing absorption, which results in a truncated domain with only trivial reflections.
+    This implementation assumes that the simulation domain without truncation is a rectangle.
+    """
+
+    def __init__(self, description: Description, round_trip: float = 1e-6, degree: int = 2):
+        self.bbox = BoundingBox2D(
+            0.0, 0.0, description.left_width + description.width + description.right_width, description.height
+        )
+        self.depth = description.absorber_depth
+        self.round_trip = round_trip
+        self.degree = degree
+        self.sigma_0 = -(self.degree + 1) * np.log(round_trip) / (2.0 * self.depth)
+
+    def eval(self, x: np.array) -> np.array:
+        """Returns modification factor to the wave number caused by the adiabatic layer.
+
+        Args:
+            x: location vector
+
+        Returns:
+            zero in areas where the adiabatic is not active, and increasingly complex values inside absorber.
+        """
+        abs_distance = self.bbox.distance(x)
+        rel_distance = abs_distance / self.depth
+
+        return self.sigma_0 * 1j * rel_distance**self.degree

src/nos/data/helmholtz/solver/wave_number_function.py

@@ -0,0 +1,11 @@
+from abc import ABC, abstractmethod
+
+import numpy as np
+
+
+class WaveNumberFunction(ABC):
+    """Abstract function that modifies the wave number on a given domain"""
+
+    @abstractmethod
+    def eval(self, x: np.array) -> np.array:
+        pass

test/data/helmholtz/solver/test_adiabatic_layer.py

@@ -0,0 +1,63 @@
+import numpy as np
+import pytest
+
+from nos.data.helmholtz.domain_properties import Description, NoneDescription
+from nos.data.helmholtz.solver import AdiabaticAbsorber
+
+
+@pytest.fixture
+def all_sides_absorber():
+    """Returns an absorber with absorber-free dimensions lower left = (0, 0), upper right = (3, 2).
+
+    Returns:
+        adiabatic absorber
+    """
+    des = Description(
+        frequencies=np.array([1]),
+        rho=1.2,
+        c=1.0,
+        left_space=0.5,
+        right_space=0.5,
+        elements=42,
+        depth=1,
+        round_trip=10,
+        directions={
+            "top": True,
+            "left": True,
+            "bottom": True,
+            "right": True,
+        },
+        crystal_description=NoneDescription("None", grid_size=1, n_x=2, n_y=2),
+    )
+    absorber = AdiabaticAbsorber(des, des.round_trip)
+
+    assert absorber.bbox.x_max == 3.0
+    assert absorber.bbox.y_max == 2.0
+
+    return absorber
+
+
+def test_inside_zero(all_sides_absorber):
+    x = np.random.random((1000, 2)) @ np.array([[3.0, 0.0], [0.0, 2.0]])  # scale to match all inside
+    assert np.allclose(all_sides_absorber.eval(x), 0)
+
+
+def test_outside_not_zero(all_sides_absorber):
+    # all outside
+    x = np.stack(
+        [
+            np.concatenate([np.random.random((500,)) - 1.1, np.random.random((500,)) + 3.1]),
+            np.concatenate([np.random.random((500,)) - 1.1, np.random.random((500,)) + 2.1]),
+        ],
+        axis=1,
+    )
+    assert all(all_sides_absorber.eval(x) != 0)
+
+
+def test_correct_value(all_sides_absorber):
+    x = np.array([[4.0, 1.0], [-1, -1]])
+
+    sigma_0 = -3 * np.log(10) / 2.0
+    sol = sigma_0 * 1j * np.array([1.0, np.sqrt(2)]) ** 2
+
+    assert np.allclose(all_sides_absorber.eval(x), sol)