Merge pull request #134 from NDF-Poli-USP/issue_0132_CE_acoustic

Add Clayton-Engquist boundary condition for acoustic wave

spyro/io/model_parameters.py

@@ -383,6 +383,13 @@ def _sanitize_absorving_boundary_condition(self):
         self.abc_pad_length = BL_obj.abc_pad_length
         self.abc_boundary_layer_type = BL_obj.abc_boundary_layer_type
 
+        self.absorb_top = dictionary.get("absorb_top", False)
+        self.absorb_bottom = dictionary.get("absorb_bottom", True)
+        self.absorb_right = dictionary.get("absorb_right", True)
+        self.absorb_left = dictionary.get("absorb_left", True)
+        self.absorb_front = dictionary.get("absorb_front", True)
+        self.absorb_back = dictionary.get("absorb_back", True)
+
     def _sanitize_output(self):
         #         default_dictionary["visualization"] = {
         #     "forward_output" : True,

spyro/solvers/acoustic_solver_construction_no_pml.py

@@ -1,5 +1,5 @@
 import firedrake as fire
-from firedrake import dx, Constant, dot, grad
+from firedrake import ds, dx, Constant, dot, grad
 
 
 def construct_solver_or_matrix_no_pml(Wave_object):
@@ -40,10 +40,27 @@ def construct_solver_or_matrix_no_pml(Wave_object):
     le = 0
     q = Wave_object.source_expression
     if q is not None:
-        le = q * v * dx(scheme=quad_rule)
+        le += q * v * dx(scheme=quad_rule)
 
     B = fire.Cofunction(V.dual())
 
+    if Wave_object.abc_active:
+        f_abc = - (1/Wave_object.c) * dot((u_n - u_nm1) / Constant(dt), v)
+        qr_s = Wave_object.surface_quadrature_rule
+        if Wave_object.absorb_top:
+            le += f_abc*ds(1, scheme=qr_s)
+        if Wave_object.absorb_bottom:
+            le += f_abc*ds(2, scheme=qr_s)
+        if Wave_object.absorb_right:
+            le += f_abc*ds(3, scheme=qr_s)
+        if Wave_object.absorb_left:
+            le += f_abc*ds(4, scheme=qr_s)
+        if Wave_object.dimension == 3:
+            if Wave_object.absorb_front:
+                le += f_abc*ds(5, scheme=qr_s)
+            if Wave_object.absorb_back:
+                le += f_abc*ds(6, scheme=qr_s)
+
     form = m1 + a - le
     lhs = fire.lhs(form)
     rhs = fire.rhs(form)

spyro/solvers/acoustic_wave.py

@@ -18,9 +18,17 @@
 )
 from ..domains.space import FE_method
 from ..utils.typing import override
+from .functionals import acoustic_energy
 
 
 class AcousticWave(Wave):
+    def __init__(self, dictionary, comm=None):
+        super().__init__(dictionary, comm=comm)
+
+        self.acoustic_energy = None
+        self.field_logger.add_functional("acoustic_energy",
+                                         lambda: fire.assemble(self.acoustic_energy))
+
     def save_current_velocity_model(self, file_name=None):
         if self.c is None:
             raise ValueError("C not loaded")
@@ -61,6 +69,8 @@ def matrix_building(self):
             self.X_np1 = fire.Function(V * Z)
             construct_solver_or_matrix_with_pml(self)
 
+        self.acoustic_energy = acoustic_energy(self)
+
     @ensemble_gradient
     def gradient_solve(self, guess=None, misfit=None, forward_solution=None):
         """Solves the adjoint problem to calculate de gradient.

spyro/solvers/functionals.py

@@ -0,0 +1,12 @@
+from firedrake import dx
+
+
+def acoustic_energy(wave):
+    '''
+    Calculates the acoustic energy as either the potential
+    energy (if u is the pressure) or the kinetic energy (if
+    u is the velocity).
+    '''
+    u = wave.get_function()
+    c = wave.c
+    return (0.5*(u/c)**2)*dx

test/test_forward_examples.py

@@ -50,6 +50,27 @@ def test_rectangle_forward():
     assert all([test1, test2, test3])
 
 
+def test_acoustic_local_abc():
+    dictionary = {}
+    dictionary["absorving_boundary_conditions"] = {
+        "status": True,
+        "damping_type": "local",
+        "absorb_top": True,
+        "absorb_bottom": True,
+        "absorb_right": True,
+        "absorb_left": True,
+    }
+    dictionary["visualization"] = {
+        "forward_output": False,
+        "acoustic_energy": True,
+        "acoustic_energy_filename": "results/acoustic_potential_energy.npy",
+    }
+    wave = spyro.examples.Camembert_acoustic(dictionary=dictionary)
+    wave.forward_solve()
+    last_acoustic_energy = wave.field_logger.get("acoustic_energy")
+    assert last_acoustic_energy < 7e-7  # The expected value was found empirically
+
+
 def test_camembert_elastic():
     from spyro.examples.camembert_elastic import wave
     wave.forward_solve()