diff --git a/src/nos/data/helmholtz/solver/__init__.py b/src/nos/data/helmholtz/solver/__init__.py index e10139ec..d14d870f 100644 --- a/src/nos/data/helmholtz/solver/__init__.py +++ b/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", ] diff --git a/src/nos/data/helmholtz/solver/adiabatic_absorber.py b/src/nos/data/helmholtz/solver/adiabatic_absorber.py new file mode 100644 index 00000000..ab076423 --- /dev/null +++ b/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 diff --git a/src/nos/data/helmholtz/solver/wave_number_function.py b/src/nos/data/helmholtz/solver/wave_number_function.py new file mode 100644 index 00000000..bebb4174 --- /dev/null +++ b/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 diff --git a/test/data/helmholtz/solver/test_adiabatic_layer.py b/test/data/helmholtz/solver/test_adiabatic_layer.py new file mode 100644 index 00000000..8c56e352 --- /dev/null +++ b/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)