started adding elastic support

cpw_elastic.py

@@ -0,0 +1,47 @@
+import numpy as np
+import spyro
+
+
+def test_cpw_calc():
+    grid_point_calculator_parameters = {
+        # Experiment parameters
+        # Here we define the frequency of the Ricker wavelet source
+        "source_frequency": 5.0,
+        "velocity_profile_type": "homogeneous",
+        "FEM_method_to_evaluate": "mass_lumped_triangle",
+        # line defines a line of point receivers with pre-established near and far
+        # offsets.
+        # Line search parameters
+        "reference_solution_file": "test/inputfiles/reference_solution_cpw.npy",
+        "reference_degree": 4,  # Degree to use in the reference case (int)
+        # grid point density to use in the reference case (float)
+        "C_reference": 5.0,
+        "desired_degree": 4,  # degree we are calculating G for. (int)
+        "C_initial": 2.2,  # Initial G for line search (float)
+        "equation_type": "isotropic_elastic",
+    }
+
+    Cpw_calc = spyro.tools.Meshing_parameter_calculator(
+        grid_point_calculator_parameters
+    )
+
+    # Check correct offset
+    source_location = Cpw_calc.initial_guess_object.source_locations[0]
+    receiver_location = Cpw_calc.initial_guess_object.receiver_locations[1]
+    sz, sx = source_location
+    rz, rx = receiver_location
+    offset = np.sqrt((sz - rz) ** 2 + (sx - rx) ** 2)
+    expected_offset_value = 2.6580067720004026
+    test1 = np.isclose(offset, expected_offset_value)
+    print(f"Checked if offset calculation is correct: {test1}")
+
+    # Check if cpw is within error TOL, starting search at min
+    min = Cpw_calc.find_minimum()
+    print(f"Minimum of {min}")
+
+    print("END")
+    # assert all([test1, test2, test3])
+
+
+if __name__ == "__main__":
+    test_cpw_calc()

spyro/tools/cells_per_wavelength_calculator.py

@@ -5,6 +5,13 @@
 import spyro
 
 
+def wave_propagator(dictionary, equation_type):
+    if equation_type == "acoustic":
+        return spyro.AcousticWave(dictionary)
+    elif equation_type == "isotropic_elastic":
+        return spyro.IsotropicWave(dictionary)
+
+
 class Meshing_parameter_calculator:
     """
     A class used to calculate meshing parameter C (cells-per-wavelength).
@@ -100,6 +107,22 @@ def __init__(self, parameters_dictionary):
             - "accepted_error_threshold": float, the accepted error threshold for the calculation
             - "desired_degree": int, the desired degree for the calculation
         """
+        parameters_dictionary.setdefault("source_frequency", 5.0)
+        parameters_dictionary.setdefault("minimum_velocity_in_the_domain", 1.5)
+        parameters_dictionary.setdefault("velocity_profile_type", "homogeneous")
+        parameters_dictionary.setdefault("velocity_model_file_name", None)
+        parameters_dictionary.setdefault("FEM_method_to_evaluate", "mass_lumped_triangle")
+        parameters_dictionary.setdefault("dimension", 2)
+        parameters_dictionary.setdefault("receiver_setup", "near")
+        parameters_dictionary.setdefault("load_reference", False)
+        parameters_dictionary.setdefault("save_reference", True)
+        parameters_dictionary.setdefault("time-step_calculation", "estimate")
+        parameters_dictionary.setdefault("C_accuracy", 0.1)
+        parameters_dictionary.setdefault("equation_type", "acoustic")
+        parameters_dictionary.setdefault("accepted_error_threshold", 0.05)
+        parameters_dictionary.setdefault("testing", False)
+
+        self.equation_type = parameters_dictionary["equation_type"]
         self.parameters_dictionary = parameters_dictionary
         self.source_frequency = parameters_dictionary["source_frequency"]
         self.minimum_velocity = parameters_dictionary[
@@ -131,10 +154,11 @@ def __init__(self, parameters_dictionary):
         self.cpw_accuracy = parameters_dictionary["C_accuracy"]
 
         # Debugging and testing  parameters
-        self._setting_up_testing_options()
+        self.reduced_obj_for_testing = self.parameters_dictionary["testing"]
 
         # Setting up reference file read or load
-        self._setting_up_reference_file_read_or_load()
+        self.save_reference = self.parameters_dictionary["save_reference"]
+        self.load_reference = self.parameters_dictionary["load_reference"]
 
         # Setting up time-step attributes
         self._setting_up_time_step()
@@ -143,26 +167,6 @@ def __init__(self, parameters_dictionary):
         self.comm = self.initial_guess_object.comm
         self.reference_solution = self.get_reference_solution()
 
-    def _setting_up_testing_options(self):
-        """
-        Sets up the testing options.
-        """
-        if "testing" in self.parameters_dictionary:
-            self.reduced_obj_for_testing = self.parameters_dictionary["testing"]
-        else:
-            self.reduced_obj_for_testing = False
-
-    def _setting_up_reference_file_read_or_load(self):
-        if "save_reference" in self.parameters_dictionary:
-            self.save_reference = self.parameters_dictionary["save_reference"]
-        else:
-            self.save_reference = False
-
-        if "load_reference" in self.parameters_dictionary:
-            self.load_reference = self.parameters_dictionary["load_reference"]
-        else:
-            self.load_reference = False
-
     def _setting_up_time_step(self):
         if "time-step_calculation" in self.parameters_dictionary:
             self.timestep_calculation = self.parameters_dictionary["time-step_calculation"]
@@ -221,7 +225,7 @@ def build_initial_guess_model(self):
         """
         dictionary = create_initial_model_for_meshing_parameter(self)
         self.initial_dictionary = dictionary
-        return spyro.AcousticWave(dictionary)
+        return wave_propagator(dictionary, self.equation_type)
 
     def get_reference_solution(self):
         """
@@ -238,9 +242,9 @@ def get_reference_solution(self):
             else:
                 filename = "reference_solution.npy"
             return np.load(filename)
-        elif self.velocity_profile_type == "heterogeneous":
+        elif self.velocity_profile_type == "heterogeneous" or self.equation_type == "isotropic_elastic":
             return self.calculate_reference_solution()
-        elif self.velocity_profile_type == "homogeneous":
+        elif self.velocity_profile_type == "homogeneous" and self.equation_type == "acoustic":
             return self.calculate_analytical_solution()
 
     def calculate_reference_solution(self):
@@ -404,7 +408,7 @@ def build_current_object(self, cpw, degree=None):
         dictionary["mesh"]["cells_per_wavelength"] = cpw
         if degree is not None:
             dictionary["options"]["degree"] = degree
-        Wave_obj = spyro.AcousticWave(dictionary)
+        Wave_obj = wave_propagator(dictionary, self.equation_type)
         if self.velocity_profile_type == "homogeneous":
             lba = self.minimum_velocity / self.source_frequency
             edge_length = lba / cpw

spyro/tools/input_models.py

@@ -131,11 +131,24 @@ def create_initial_model_for_meshing_parameter(Meshing_calc_obj):
     """
     dimension = Meshing_calc_obj.dimension
     if dimension == 2:
-        return create_initial_model_for_meshing_parameter_2D(Meshing_calc_obj)
+        model_dictionary = create_initial_model_for_meshing_parameter_2D(Meshing_calc_obj)
     elif dimension == 3:
-        return create_initial_model_for_meshing_parameter_3D(Meshing_calc_obj)
+        model_dictionary = create_initial_model_for_meshing_parameter_3D(Meshing_calc_obj)
     else:
         raise ValueError("Dimension is not 2 or 3")
+
+    if Meshing_calc_obj.equation_type == "isotropic_elastic" and dimension == 2:
+        model_dictionary = add_elastic_to_dictionary(model_dictionary)
+
+    return model_dictionary
+
+
+def add_elastic_to_dictionary(dictionary):
+    model_dictionary = dictionary
+    dictionary["absorving_boundary_conditions"]["damping_type"] = "local"
+    dictionary["absorving_boundary_conditions"]["pad"] = 0.0
+
+    return model_dictionary
 
 
 def create_initial_model_for_meshing_parameter_2D(Meshing_calc_obj):