diff --git a/spyro/io/model_parameters.py b/spyro/io/model_parameters.py index df368466..146be60f 100644 --- a/spyro/io/model_parameters.py +++ b/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, diff --git a/spyro/solvers/acoustic_solver_construction_no_pml.py b/spyro/solvers/acoustic_solver_construction_no_pml.py index 7ce6a9ad..36a6e43b 100644 --- a/spyro/solvers/acoustic_solver_construction_no_pml.py +++ b/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) diff --git a/spyro/solvers/acoustic_wave.py b/spyro/solvers/acoustic_wave.py index 2f13b157..8406e34a 100644 --- a/spyro/solvers/acoustic_wave.py +++ b/spyro/solvers/acoustic_wave.py @@ -17,6 +17,7 @@ ) from ..domains.space import FE_method from ..utils.typing import override +from .functionals import acoustic_energy try: from SeismicMesh import write_velocity_model @@ -26,6 +27,13 @@ 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") @@ -66,6 +74,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. diff --git a/spyro/solvers/functionals.py b/spyro/solvers/functionals.py new file mode 100644 index 00000000..16c1f172 --- /dev/null +++ b/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 diff --git a/test/test_forward_examples.py b/test/test_forward_examples.py index 77e89062..deb9394f 100644 --- a/test/test_forward_examples.py +++ b/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()