Merge pull request #126 from Olender/issue_0041

New features of isotropic elastic wave propagator

.gitignore

@@ -7,3 +7,5 @@
 *.hdf5
 *.msh
 *.segy
+*.png
+*.vtk

docker/Dockerfile

@@ -1,4 +1,4 @@
-FROM firedrakeproject/firedrake:2023-09 AS spyro_base
+FROM firedrakeproject/firedrake:2024-09 AS spyro_base
 
 USER root
 RUN apt-get update \
@@ -9,8 +9,6 @@ RUN apt-get update \
 USER firedrake
 RUN . /home/firedrake/firedrake/bin/activate; pip3 install wheel --upgrade
 RUN . /home/firedrake/firedrake/bin/activate; pip3 install scipy
-RUN . /home/firedrake/firedrake/bin/activate; pip3 install roltrilinos
-RUN . /home/firedrake/firedrake/bin/activate; pip3 install ROL
 RUN . /home/firedrake/firedrake/bin/activate; pip3 install matplotlib
 RUN . /home/firedrake/firedrake/bin/activate; pip3 install segyio
 
@@ -19,16 +17,21 @@ USER root
 RUN apt-get update && apt-get install -y libgmp3-dev libmpfr-dev libcgal-dev python3-tk
 USER firedrake
 # The PyPI SeismicMesh has bugs. It is better to install from GitHub
-RUN . /home/firedrake/firedrake/bin/activate; pip3 install --no-dependencies git+https://github.com/krober10nd/SeismicMesh.git
+RUN . /home/firedrake/firedrake/bin/activate; pip3 install --no-dependencies git+https://github.com/NDF-Poli-USP/SeismicMesh.git
 RUN . /home/firedrake/firedrake/bin/activate; pip3 install pyamg
 RUN . /home/firedrake/firedrake/bin/activate; pip3 install meshio
+RUN . /home/firedrake/firedrake/bin/activate; pip3 install siphash24
 
 FROM spyro_base AS spyro_release
 
 RUN . /home/firedrake/firedrake/bin/activate; pip install git+https://github.com/Olender/spyro-1.git
 
 FROM spyro_base AS spyro_development
 
+# For notebook development
+RUN . /home/firedrake/firedrake/bin/activate; pip3 install notebook
+EXPOSE 8888
+
 RUN . /home/firedrake/firedrake/bin/activate; pip3 install pytest
 WORKDIR /home/firedrake
-RUN echo "/home/firedrake/shared/spyro" >> /home/firedrake/firedrake/lib/python3.10/site-packages/shared.pth
+RUN echo "/home/firedrake/shared/spyro" >> /home/firedrake/firedrake/lib/python3.12/site-packages/shared.pth

docker/README.md

@@ -36,5 +36,7 @@ docker run -v $PWD:/home/firedrake/shared/spyro -it devtag:1.0
 For running the automated tests:
 ````
 cd shared/spyro/
-python3 -m pytest --maxfail=1 test/
+python3 -m pytest test/
+python3 -m pytest test_3d/
+mpiexec -n 6 python3 -m pytest test_parallel/
 ````

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.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/__init__.py

@@ -5,6 +5,7 @@
 from .sources.Sources import Sources, ricker_wavelet, full_ricker_wavelet
 from .solvers.wave import Wave
 from .solvers.acoustic_wave import AcousticWave
+from .solvers.elastic_wave.isotropic_wave import IsotropicWave
 from .solvers.inversion import FullWaveformInversion
 
 # from .solvers.dg_wave import DG_Wave
@@ -50,4 +51,5 @@
     "RectangleMesh",
     "PeriodicRectangleMesh",
     "BoxMesh",
+    "IsotropicWave",
 ]

spyro/domains/space.py

@@ -1,7 +1,7 @@
-from firedrake import *  # noqa:F403
+from firedrake import (FiniteElement, FunctionSpace, VectorElement)
 
 
-def FE_method(mesh, method, degree):
+def FE_method(mesh, method, degree, dim=1):
     """Define the finite element space:
 
     Parameters:
@@ -12,6 +12,8 @@ def FE_method(mesh, method, degree):
         Method to be used for the finite element space.
     degree: int
         Degree of the finite element space.
+    dim: int
+        Number of degrees of freedom per node.
 
     Returns:
     --------
@@ -20,21 +22,24 @@ def FE_method(mesh, method, degree):
     """
 
     if method == "mass_lumped_triangle":
-        element = FiniteElement(  # noqa: F405
+        element = FiniteElement(
             "KMV", mesh.ufl_cell(), degree=degree, variant="KMV"
         )
     elif method == "spectral_quadrilateral":
-        element = FiniteElement(  # noqa: F405
+        element = FiniteElement(
             "CG", mesh.ufl_cell(), degree=degree, variant="spectral"
         )
     elif method == "DG_triangle" or "DG_quadrilateral" or "DG":
         element = FiniteElement(
             "DG", mesh.ufl_cell(), degree=degree
-        )  # noqa: F405
+        )
     elif method == "CG_triangle" or "CG_quadrilateral" or "CG":
         element = FiniteElement(
             "CG", mesh.ufl_cell(), degree=degree
-        )  # noqa: F405
+        )
 
-    function_space = FunctionSpace(mesh, element)  # noqa: F405
+    if dim > 1:
+        element = VectorElement(element, dim=dim)
+
+    function_space = FunctionSpace(mesh, element)
     return function_space

spyro/examples/camembert_elastic.py

@@ -0,0 +1,99 @@
+'''
+From Mohammad Mehdi Ghorbani's dissertation:
+https://www.teses.usp.br/teses/disponiveis/3/3152/tde-16032023-085210/pt-br.php
+'''
+import firedrake as fire
+import numpy as np
+import spyro
+
+L = 500 # [m]
+
+rho = 7850          # [kg/m3]
+lambda_in = 6.86e9  # [Pa]
+lambda_out = 9.88e9 # [Pa]
+mu_in = 3.86e9      # [Pa]
+mu_out = 5.86e9     # [Pa]
+
+smag = 1e6
+freq = 2 # Central frequency of Ricker wavelet [Hz]
+hf =  90 # [m]
+hs = 100 # [m]
+source_locations = spyro.create_transect((-hf, 0.2*L), (-hf, 0.8*L), 3)
+receiver_locations = spyro.create_transect((-hs, 0), (-hs, L), 40)
+source_locations = [[-hf, 0.5*L]]
+
+time_step = 2e-4 # [s]
+final_time = 1.5 # [s]
+out_freq = int(0.01/time_step)
+
+n = 20
+mesh = fire.RectangleMesh(n, n, 0, L, originX=-L, diagonal='crossed')
+z, x = fire.SpatialCoordinate(mesh)
+
+zc = 250 # [m]
+xc = 250 # [m]
+ri =  50 # [m]
+camembert = lambda v_inside, v_outside: fire.conditional(
+    (z - zc) ** 2 + (x - xc) ** 2 < ri**2, v_inside, v_outside)
+
+d = {}
+
+d["options"] = {
+    "cell_type": "T",
+    "variant": "lumped",
+    "degree": 4,
+    "dimension": 2,
+}
+
+d["parallelism"] = {
+    "type": "automatic",
+}
+
+d["mesh"] = {
+    "user_mesh": mesh,
+}
+
+d["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": source_locations,
+    "frequency": freq,
+    "delay": 0,
+    "delay_type": "time",
+    "amplitude": np.array([0, smag]),
+    #"amplitude": smag * np.eye(2),
+    #"amplitude": smag * np.array([[0, 1], [-1, 0]]),
+    "receiver_locations": receiver_locations,
+}
+
+d["synthetic_data"] = {
+    "type": "object",
+    "density": fire.Constant(rho),
+    "lambda": camembert(lambda_in, lambda_out),
+    "mu": camembert(mu_in, mu_out),
+    "real_velocity_file": None,
+}
+
+d["time_axis"] = {
+    "initial_time": 0,
+    "final_time": final_time,
+    "dt": time_step,
+    "output_frequency": out_freq,
+    "gradient_sampling_frequency": 1,
+}
+
+d["visualization"] = {
+    "forward_output": True,
+    "forward_output_filename": "results/forward_output.pvd",
+    "fwi_velocity_model_output": False,
+    "gradient_output": False,
+    "adjoint_output": False,
+    "debug_output": False,
+}
+
+d["absorving_boundary_conditions"] = {
+    "status": True,
+    "damping_type": "local",
+}
+
+wave = spyro.IsotropicWave(d)
+wave.set_mesh(user_mesh=mesh, mesh_parameters={})

spyro/examples/elastic_cube_3D.py

@@ -0,0 +1,90 @@
+'''
+From Turkel et al. (2023):
+https://doi.org/10.1016/j.wavemoti.2022.103109
+'''
+import numpy as np
+import spyro
+
+L = 2000   # Size of the edges of the cube [m]
+N = 5      # Number of elements in each direction
+h = L/N    # Element size [m]
+
+c_p = 5000 # P-wave velocity [m/s]
+c_s = 2500 # S-wave velocity [m/s]
+rho = 1000 # Density [kg/m3]
+
+smag = 1e9 # Source magnitude
+freq = 1   # Source frequency [Hz]
+
+final_time = 2
+time_step = 5e-4
+out_freq = (0.01/time_step)
+
+print(f'Element size: {h} m')
+print(f'Cross time  : {h/c_p} s')
+print(f'Time step   : {time_step} s')
+
+d = {}
+
+d["options"] = {
+    "cell_type": "T",
+    "variant": "lumped",
+    "degree": 3,
+    "dimension": 3,
+}
+
+d["parallelism"] = {
+    "type": "automatic",
+}
+
+d["mesh"] = {
+    "Lz": L,
+    "Lx": L,
+    "Ly": L,
+    "h": h,
+    "mesh_file": None,
+    "mesh_type": "firedrake_mesh",
+}
+
+d["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": [[-L/2, L/2, L/2]],
+    "frequency": freq,
+    "delay": 0,
+    "delay_type": "time",
+    "amplitude": np.array([smag, 0, 0]),
+    "receiver_locations": [],
+}
+
+d["synthetic_data"] = {
+    "type": "object",
+    "density": rho,
+    "p_wave_velocity": c_p,
+    "s_wave_velocity": c_s,
+    "real_velocity_file": None,
+}
+
+d["time_axis"] = {
+    "initial_time": 0,
+    "final_time": final_time,
+    "dt": time_step,
+    "output_frequency": out_freq,
+    "gradient_sampling_frequency": 1,
+}
+
+d["visualization"] = {
+    "forward_output": True,
+    "forward_output_filename": "results/elastic_cube_3D.pvd",
+    "fwi_velocity_model_output": False,
+    "gradient_output": False,
+    "adjoint_output": False,
+    "debug_output": False,
+}
+
+d["absorving_boundary_conditions"] = {
+    "status": True,
+    "damping_type": "local",
+}
+
+wave = spyro.IsotropicWave(d)
+wave.set_mesh(mesh_parameters={'dx': h})

spyro/examples/rectangle.py

@@ -181,26 +181,3 @@ def multiple_layer_velocity_model(self, z_switch, layers):
                 )
         # cond = fire.conditional(self.mesh_z > z_switch, layer1, layer2)
         self.set_initial_velocity_model(conditional=cond)
-
-
-# class Rectangle(AcousticWave):
-#     def __init__(self, model_dictionary=None, comm=None):
-#         model_parameters = Rectangle_parameters(
-#             dictionary=model_dictionary, comm=comm
-#         )
-#         super().__init__(
-#             model_parameters=model_parameters, comm=model_parameters.comm
-#         )
-#         comm = self.comm
-#         num_sources = self.number_of_sources
-#         if comm.comm.rank == 0 and comm.ensemble_comm.rank == 0:
-#             print(
-#                 "INFO: Distributing %d shot(s) across %d core(s). \
-#                     Each shot is using %d cores"
-#                 % (
-#                     num_sources,
-#                     fire.COMM_WORLD.size,
-#                     fire.COMM_WORLD.size / comm.ensemble_comm.size,
-#                 ),
-#                 flush=True,
-#             )

spyro/io/boundary_layer_io.py

@@ -45,7 +45,7 @@ class read_boundary_layer:
 
     def __init__(self, abc_dictionary):
         self.dictionary = abc_dictionary
-        if self.dictionary["status"] is False:
+        if self.dictionary["status"] is False or self.dictionary["damping_type"] == "local":
             self.abc_exponent = None
             self.abc_cmax = None
             self.abc_R = None