Issue 0041 - Add elastic wave propagation

notebook_tutorials/altering_time_integration.ipynb

@@ -338,7 +338,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Wave_obj._get_initial_velocity_model()\n",
+    "Wave_obj._initialize_model_parameters()\n",
     "Wave_obj.c = Wave_obj.initial_velocity_model\n",
     "Wave_obj.matrix_building()"
    ]
@@ -347,7 +347,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The `_get_initial_velocity_model` is not actually necessary for this specific case (where the velocity model is already created in the object) and does not do anything in this case. However, it needs to be called, because in larger cases the velocity model is only loaded, by this method, into our function space right before it is needed (to conserve memory). Spyro will build the solver object based on the velocity model that the `c` attribute points to. This can be diferent than the `initial_velocity_model`, especially in inversion problems."
+    "The `_initialize_model_parameters` is not actually necessary for this specific case (where the velocity model is already created in the object) and does not do anything in this case. However, it needs to be called, because in larger cases the velocity model is only loaded, by this method, into our function space right before it is needed (to conserve memory). Spyro will build the solver object based on the velocity model that the `c` attribute points to. This can be diferent than the `initial_velocity_model`, especially in inversion problems."
    ]
   },
   {

notebook_tutorials/simple_forward_with_overthrust.html

@@ -8051,7 +8051,7 @@ <h1 id="Simple-forward-with-overthrust">Simple forward with overthrust<a class="
 <div class="jp-InputArea jp-Cell-inputArea"><div class="jp-InputPrompt jp-InputArea-prompt">
 </div><div class="jp-RenderedHTMLCommon jp-RenderedMarkdown jp-MarkdownOutput" data-mime-type="text/markdown">
 <p>From the plot above, we can notice two things. First, it is rotated by 90 degrees. In cases like this, the saved PNG image will be in the correct orientation. However, the plotted image in this notebook is rotated because the velocity model we used, commonly with segy files, has the Z-axis before the X-axis in data input.</p>
-<p>Another observation is that, unlike in the <strong>simple forward tutorial</strong>, we looked at the experiment layout after running the forward solve. For memory purposes, a 2D or 3D velocity file is only actually interpolated into our domain when it is necessary for another method of the Wave object class. If you need to force the interpolation sooner, call the _get_initial_velocity_model() method.</p>
+<p>Another observation is that, unlike in the <strong>simple forward tutorial</strong>, we looked at the experiment layout after running the forward solve. For memory purposes, a 2D or 3D velocity file is only actually interpolated into our domain when it is necessary for another method of the Wave object class. If you need to force the interpolation sooner, call the _initialize_model_parameters() method.</p>
 <p>It is also important to note that even though receivers look like a line, they are actually located in points, which can be visible by zooming into the image, not coinciding with nodes.</p>
 </div>
 </div>

notebook_tutorials/simple_forward_with_overthrust.ipynb

@@ -397,7 +397,7 @@
    "source": [
     "From the plot above, we can notice two things. First, it is rotated by 90 degrees. In cases like this, the saved PNG image will be in the correct orientation. However, the plotted image in this notebook is rotated because the velocity model we used, commonly with segy files, has the Z-axis before the X-axis in data input.\n",
     "\n",
-    "Another observation is that, unlike in the **simple forward tutorial**, we looked at the experiment layout after running the forward solve. For memory purposes, a 2D or 3D velocity file is only actually interpolated into our domain when it is necessary for another method of the Wave object class. If you need to force the interpolation sooner, call the _get_initial_velocity_model() method.\n",
+    "Another observation is that, unlike in the **simple forward tutorial**, we looked at the experiment layout after running the forward solve. For memory purposes, a 2D or 3D velocity file is only actually interpolated into our domain when it is necessary for another method of the Wave object class. If you need to force the interpolation sooner, call the _initialize_model_parameters() method.\n",
     "\n",
     "It is also important to note that even though receivers look like a line, they are actually located in points, which can be visible by zooming into the image, not coinciding with nodes."
    ]

spyro/io/model_parameters.py

@@ -608,9 +608,9 @@ def _sanitize_optimization_and_velocity(self):
             if "velocity_conditional" not in dictionary["synthetic_data"]:
                 self.velocity_model_type = None
                 warnings.warn(
-                    "No velocity model set initially. If using \
-                        user defined conditional or expression, please \
-                            input it in the Wave object."
+                    "No velocity model set initially. If using " \
+                        "user defined conditional or expression, please " \
+                            "input it in the Wave object."
                 )
 
         if "velocity_conditional" in dictionary["synthetic_data"]:

spyro/solvers/acoustic_wave.py

@@ -40,7 +40,7 @@ def forward_solve(self):
         if self.function_space is None:
             self.force_rebuild_function_space()
 
-        self._get_initial_velocity_model()
+        self._initialize_model_parameters()
         self.c = self.initial_velocity_model
         self.matrix_building()
         self.wave_propagator()
@@ -144,3 +144,33 @@ def reset_pressure(self):
             self.u_n.assign(0.0)
         except:
             warnings.warn("No pressure to reset")
+
+    #@override
+    def _initialize_model_parameters(self):
+        if self.initial_velocity_model is not None:
+            return None
+
+        if self.initial_velocity_model_file is None:
+            raise ValueError("No velocity model or velocity file to load.")
+
+        if self.initial_velocity_model_file.endswith(".segy"):
+            vp_filename, vp_filetype = os.path.splitext(
+                self.initial_velocity_model_file
+            )
+            warnings.warn("Converting segy file to hdf5")
+            write_velocity_model(
+                self.initial_velocity_model_file, ofname=vp_filename
+            )
+            self.initial_velocity_model_file = vp_filename + ".hdf5"
+
+        if self.initial_velocity_model_file.endswith((".hdf5", ".h5")):
+            self.initial_velocity_model = interpolate(
+                self,
+                self.initial_velocity_model_file,
+                self.function_space.sub(0),
+            )
+
+        if self.debug_output:
+            fire.File("initial_velocity_model.pvd").write(
+                self.initial_velocity_model, name="velocity"
+            )

spyro/solvers/elastic_wave/elastic_wave.py

@@ -0,0 +1,29 @@
+from abc import abstractmethod
+
+from ..wave import Wave
+
+class ElasticWave(Wave):
+    '''Base class for elastic wave propagators'''
+    def __init__(self, dictionary, comm=None):
+        super().__init__(dictionary, comm=comm)
+
+    #@override
+    def _initialize_model_parameters(self):
+        d = self.input_dictionary.get("synthetic_data", False)
+        if bool(d) and "type" in d:
+            if d["type"] == "object":
+                self.initialize_model_parameters_from_object(d)
+            elif d["type"] == "file":
+                self.initialize_model_parameters_from_file(d)
+            else:
+                raise Exception(f"Invalid synthetic data type: {d['type']}")
+        else:
+            raise Exception("Input dictionary must contain ['synthetic_data']['type']")
+
+    @abstractmethod
+    def initialize_model_parameters_from_object(self, synthetic_data_dict):
+        pass
+
+    @abstractmethod
+    def initialize_model_parameters_from_file(self, synthetic_data_dict):
+        pass

spyro/solvers/elastic_wave/isotropic_wave.py

@@ -0,0 +1,48 @@
+from .elastic_wave import ElasticWave
+
+class IsotropicWave(ElasticWave):
+    '''Isotropic elastic wave propagator'''
+    def __init__(self, dictionary, comm=None):
+        super().__init__(dictionary, comm=comm)
+
+        self.rho = None   # Density
+        self.lmbda = None # First Lame parameter
+        self.mu = None    # Second Lame parameter
+        self.c_s = None   # Secondary wave velocity
+
+    #@override
+    def initialize_model_parameters_from_object(self, synthetic_data_dict: dict):
+        self.rho = synthetic_data_dict.get("density", None)
+        self.lmbda = synthetic_data_dict.get("lambda", 
+                                            synthetic_data_dict.get("lame_first", None))
+        self.mu = synthetic_data_dict.get("mu", 
+                                          synthetic_data_dict.get("lame_second", None))
+        self.c = synthetic_data_dict.get("p_wave_velocity", None)
+        self.c_s = synthetic_data_dict.get("s_wave_velocity", None)
+
+        # Check if {rho, lambda, mu} is set and {c, c_s} are not
+        option_1 = bool(self.rho) and \
+                   bool(self.lmbda) and \
+                   bool(self.mu) and \
+                   not bool(self.c) and \
+                   not bool(self.c_s)
+        # Check if {rho, c, c_s} is set and {lambda, mu} are not
+        option_2 = bool(self.rho) and \
+                   bool(self.c) and \
+                   bool(self.c_s) and \
+                   not bool(self.lmbda) and \
+                   not bool(self.mu)
+
+        if not option_1 and not option_2:
+            raise Exception(f"Inconsistent selection of isotropic elastic wave parameters:\n" \
+                            f"    Density        : {bool(self.rho)}\n"\
+                            f"    Lame first     : {bool(self.lmbda)}\n"\
+                            f"    Lame second    : {bool(self.mu)}\n"\
+                            f"    P-wave velocity: {bool(self.c)}\n"\
+                            f"    S-wave velocity: {bool(self.c_s)}\n"\
+                            "The valid options are \{Density, Lame first, Lame second\} "\
+                            "or \{Density, P-wave velocity, S-wave velocity\}")
+
+    #@override
+    def initialize_model_parameters_from_file(self, synthetic_data_dict):
+        raise NotImplementedError

spyro/solvers/time_integration_central_difference.py

@@ -79,8 +79,8 @@ def central_difference(Wave_object, source_id=0):
         if (step - 1) % Wave_object.output_frequency == 0:
             assert (
                 fire.norm(u_n) < 1
-            ), "Numerical instability. Try reducing dt or building the \
-                mesh differently"
+            ), "Numerical instability. Try reducing dt or building the " \
+               "mesh differently"
             if Wave_object.forward_output:
                 output.write(u_n, time=t, name="Pressure")
 
@@ -184,8 +184,8 @@ def mixed_space_central_difference(Wave_object, source_id=0):
         if (step - 1) % Wave_object.output_frequency == 0:
             assert (
                 fire.norm(X_np1.sub(0)) < 1
-            ), "Numerical instability. Try reducing dt or building the \
-                mesh differently"
+            ), "Numerical instability. Try reducing dt or building the " \
+               "mesh differently"
             if Wave_object.forward_output:
                 output.write(X_np1.sub(0), time=t, name="Pressure")
 
@@ -292,8 +292,8 @@ def central_difference_MMS(Wave_object, source_id=0):
         if (step - 1) % Wave_object.output_frequency == 0:
             assert (
                 fire.norm(u_n) < 1
-            ), "Numerical instability. Try reducing dt or building the \
-                mesh differently"
+            ), "Numerical instability. Try reducing dt or building the " \
+               "mesh differently"
             if Wave_object.forward_output:
                 output.write(u_n, time=t, name="Pressure")
             if t > 0:

spyro/solvers/wave.py

@@ -209,8 +209,8 @@ def set_initial_velocity_model(
             self.initial_velocity_model = vp
         else:
             raise ValueError(
-                "Please specify either a conditional, expression, firedrake \
-                    function or new file name (segy or hdf5)."
+                "Please specify either a conditional, expression, firedrake " \
+                    "function or new file name (segy or hdf5)."
             )
         if output:
             fire.File("initial_velocity_model.pvd").write(
@@ -223,35 +223,10 @@ def _map_sources_and_receivers(self):
         else:
             self.sources = None
         self.receivers = Receivers(self)
-
-    def _get_initial_velocity_model(self):
-        if self.initial_velocity_model is not None:
-            return None
-
-        if self.initial_velocity_model_file is None:
-            raise ValueError("No velocity model or velocity file to load.")
-
-        if self.initial_velocity_model_file.endswith(".segy"):
-            vp_filename, vp_filetype = os.path.splitext(
-                self.initial_velocity_model_file
-            )
-            warnings.warn("Converting segy file to hdf5")
-            write_velocity_model(
-                self.initial_velocity_model_file, ofname=vp_filename
-            )
-            self.initial_velocity_model_file = vp_filename + ".hdf5"
-
-        if self.initial_velocity_model_file.endswith(".hdf5"):
-            self.initial_velocity_model = interpolate(
-                self,
-                self.initial_velocity_model_file,
-                self.function_space.sub(0),
-            )
-
-        if self.debug_output:
-            fire.File("initial_velocity_model.pvd").write(
-                self.initial_velocity_model, name="velocity"
-            )
+
+    @abstractmethod
+    def _initialize_model_parameters(self):
+        pass
 
     def _build_function_space(self):
         self.function_space = FE_method(self.mesh, self.method, self.degree)

spyro/tools/cells_per_wavelength_calculator.py

@@ -347,7 +347,7 @@ def find_minimum(self, starting_cpw=None, TOL=None, accuracy=None, savetxt=False
 
             # Running forward model
             Wave_obj = self.build_current_object(cpw)
-            Wave_obj._get_initial_velocity_model()
+            Wave_obj._initialize_model_parameters() # TO REVIEW: call to protected method
 
             # Setting up time-step
             if self.timestep_calculation != "float":

test/test_isotropic_wave.py

@@ -0,0 +1,67 @@
+import pytest
+
+from spyro.solvers.elastic_wave.isotropic_wave import IsotropicWave
+
+# TO REVIEW: it is extra work to have to define this dictionary everytime
+# Here I listed only the required parameters for running to get a view of
+# what is currently necessary. Note that the dictionary is not even complete
+dummy_dict = {
+    "options": {
+        "cell_type": "T",
+        "variant": "lumped",
+    },
+    "time_axis": {
+        "final_time": 1,
+        "dt": 0.001,
+        "output_frequency": 100,
+        "gradient_sampling_frequency": 1,
+    },
+    "mesh": {},
+    "acquisition": {
+        "receiver_locations": [],
+        "source_type": "ricker",
+        "source_locations": [(0, 0)],
+        "frequency": 5.0,
+    },
+}
+
+def test_initialize_model_parameters_from_object_missing_parameters():
+    synthetic_dict = {
+        "type": "object",
+    }
+    wave = IsotropicWave(dummy_dict)
+    with pytest.raises(Exception) as e:
+        wave.initialize_model_parameters_from_object(synthetic_dict)
+
+def test_initialize_model_parameters_from_object_first_option():
+    synthetic_dict = {
+        "type": "object",
+        "density": 1,
+        "lambda": 2,
+        "mu": 3,
+    }
+    wave = IsotropicWave(dummy_dict)
+    wave.initialize_model_parameters_from_object(synthetic_dict)
+
+def test_initialize_model_parameters_from_object_second_option():
+    synthetic_dict = {
+        "type": "object",
+        "density": 1,
+        "p_wave_velocity": 2,
+        "s_wave_velocity": 3,
+    }
+    wave = IsotropicWave(dummy_dict)
+    wave.initialize_model_parameters_from_object(synthetic_dict)
+
+def test_initialize_model_parameters_from_object_redundant():
+    synthetic_dict = {
+        "type": "object",
+        "density": 1,
+        "lmbda": 2,
+        "mu": 3,
+        "p_wave_velocity": 2,
+        "s_wave_velocity": 3,
+    }
+    wave = IsotropicWave(dummy_dict)
+    with pytest.raises(Exception) as e:
+        wave.initialize_model_parameters_from_object(synthetic_dict)