Merge pull request #137 from NDF-Poli-USP/issue_0133-create-new-spyro…

…-version-were-every-dependency-not-pip-installable-is-optional

Issue 0133 create new spyro version were every dependency not pip installable is optional

setup.py

@@ -2,16 +2,17 @@
 
 setup(
     name="spyro",
-    version="0.2.0",
-    license="GPL v3",
-    description="acoustic wave modeling with the finite element method",
-    author="Keith J. Roberts, Alexandre F. G. Olender, Lucas Franceschini",
-    url="https://github.com/krober10nd/spyro",
+    version="0.9.0",
+    license="LGPL v3",
+    description="Wave modeling with the finite element method",
+    author="Keith J. Roberts, Alexandre F. G. Olender, Eduardo Moscatelli de Souza, Daiane I. Dolci, Thiago Dias dos Santos, Lucas Franceschini",
+    url="https://github.com/NDF-Poli-USP/spyro",
     packages=find_packages(),
     install_requires=[
+        "firedrake",
         "numpy",
         "scipy",
         "matplotlib",
-        "exdown==0.7.0",
-        "segyio"],
+        "segyio",
+        "meshio"],
 )

spyro/meshing/meshing_functions.py

@@ -1,7 +1,11 @@
 import firedrake as fire
-import SeismicMesh
 import meshio
 
+try:
+    import SeismicMesh
+except ImportError:
+    SeismicMesh = None
+
 
 def cells_per_wavelength(method, degree, dimension):
     cell_per_wavelength_dictionary = {
@@ -13,9 +17,9 @@ def cells_per_wavelength(method, degree, dimension):
         'mlt3tet': 3.72,
     }
 
-    if dimension == 2 and (method == 'KMV' or method == 'CG'):
+    if dimension == 2 and (method == 'MLT' or method == 'CG'):
         cell_type = 'tri'
-    if dimension == 3 and (method == 'KMV' or method == 'CG'):
+    if dimension == 3 and (method == 'MLT' or method == 'CG'):
         cell_type = 'tet'
 
     key = method.lower()+str(degree)+cell_type
@@ -221,6 +225,8 @@ def set_seismicmesh_parameters(
         -------
         None
         """
+        if SeismicMesh is None:
+            raise ImportError("SeismicMesh is not available. Please install it to use this function.")
         if cpw is not None:
             self.cpw = cpw
         if velocity_model is not None:
@@ -258,6 +264,8 @@ def create_mesh(self):
         if self.mesh_type == "firedrake_mesh":
             return self.create_firedrake_mesh()
         elif self.mesh_type == "SeismicMesh":
+            if SeismicMesh is None:
+                raise ImportError("SeismicMesh is not available. Please install it to use this function.")
             return self.create_seismicmesh_mesh()
         else:
             raise ValueError("mesh_type is not supported")

spyro/solvers/acoustic_wave.py

@@ -1,7 +1,6 @@
 import firedrake as fire
 import warnings
 import os
-from SeismicMesh import write_velocity_model
 
 from .wave import Wave
 
@@ -20,6 +19,12 @@
 from ..utils.typing import override
 from .functionals import acoustic_energy
 
+try:
+    from SeismicMesh import write_velocity_model
+    SEISMIC_MESH_AVAILABLE = True
+except ImportError:
+    SEISMIC_MESH_AVAILABLE = False
+
 
 class AcousticWave(Wave):
     def __init__(self, dictionary, comm=None):
@@ -109,6 +114,8 @@ def _initialize_model_parameters(self):
                 raise ValueError("No velocity model or velocity file to load.")
 
             if self.initial_velocity_model_file.endswith(".segy"):
+                if not SEISMIC_MESH_AVAILABLE:
+                    raise ImportError("SeismicMesh is required to convert segy files.")
                 vp_filename, vp_filetype = os.path.splitext(
                     self.initial_velocity_model_file
                 )

spyro/tools/velocity_smoother.py

@@ -70,6 +70,7 @@ def smooth_velocity_field_file(input_filename, output_filename, sigma, show=Fals
         plt.xlabel("x-direction (m)")
         plt.ylabel("z-direction (m)")
         ax.axis("equal")
+        plt.savefig(output_filename+".png")
         plt.show()
 
     return None

test/test_cpw_calc.py

@@ -1,8 +1,17 @@
 import numpy as np
+import pytest
 import spyro
 
 
-def test_cpw_calc():
+def is_seismicmesh_installed():
+    try:
+        import SeismicMesh  # noqa: F401
+        return True
+    except ImportError:
+        return False
+
+
+def run_test_cpw_calc(FEM_method_to_evaluate, correct_cpw):
     grid_point_calculator_parameters = {
         # Experiment parameters
         # Here we define the frequency of the Ricker wavelet source
@@ -17,7 +26,7 @@ def test_cpw_calc():
         "velocity_model_file_name": None,
         # FEM to evaluate such as `KMV` or `spectral`
         # (GLL nodes on quads and hexas)
-        "FEM_method_to_evaluate": "mass_lumped_triangle",
+        "FEM_method_to_evaluate": FEM_method_to_evaluate,
         "dimension": 2,  # Domain dimension. Either 2 or 3.
         # Either near or line. Near defines a receiver grid near to the source,
         "receiver_setup": "near",
@@ -61,11 +70,25 @@ def test_cpw_calc():
     # Check if cpw is within error TOL, starting search at min
     min = Cpw_calc.find_minimum()
     print(f"Minimum of {min}")
-    test3 = np.isclose(2.3, min)
+    test3 = np.isclose(correct_cpw, min)
 
     print("END")
     assert all([test1, test2, test3])
 
 
+@pytest.mark.skipif(not is_seismicmesh_installed(), reason="SeismicMesh is not installed")
+def test_cpw_calc_triangles():
+    method = "mass_lumped_triangle"
+    correct_cpw = 2.3
+    return run_test_cpw_calc(method, correct_cpw)
+
+
+def test_cpw_calc_quads():
+    method = "spectral_quadrilateral"
+    correct_cpw = 2.5
+    return run_test_cpw_calc(method, correct_cpw)
+
+
 if __name__ == "__main__":
-    test_cpw_calc()
+    test_cpw_calc_triangles()
+    test_cpw_calc_quads()

test/test_elastic_local_abc.py

@@ -6,26 +6,43 @@
 expected_mechanical_energy = 0.25
 
 
+def has_sufficient_memory():
+    meminfo = {}
+    with open('/proc/meminfo') as f:
+        for line in f:
+            parts = line.split(':')
+            if len(parts) == 2:
+                meminfo[parts[0].strip()] = parts[1].strip()
+    total_memory_kb = int(meminfo.get('MemTotal', '0 kB').split()[0])
+    total_memory_gb = total_memory_kb / 1024 / 1024
+    print(f"Total system memory {total_memory_gb}")
+    return total_memory_gb > 16
+
+
+@pytest.mark.skipif(not has_sufficient_memory(), reason="Insufficient memory")
 def test_stacey_abc():
     wave = build_solver("Stacey", "backward")
     wave.forward_solve()
     last_mechanical_energy = wave.field_logger.get("mechanical_energy")
     assert last_mechanical_energy < expected_mechanical_energy
 
 
+@pytest.mark.skipif(not has_sufficient_memory(), reason="Insufficient memory")
 def test_clayton_engquist_abc():
     wave = build_solver("CE_A1", "backward")
     wave.forward_solve()
     last_mechanical_energy = wave.field_logger.get("mechanical_energy")
     assert last_mechanical_energy < expected_mechanical_energy
 
 
+@pytest.mark.skipif(not has_sufficient_memory(), reason="Insufficient memory")
 def test_with_central():
     wave = build_solver("Stacey", "central")
     with pytest.raises(AssertionError):
         wave.forward_solve()
 
 
+@pytest.mark.skipif(not has_sufficient_memory(), reason="Insufficient memory")
 def test_with_backward_2nd():
     wave = build_solver("Stacey", "backward_2nd")
     wave.forward_solve()

test/test_meshing.py

@@ -0,0 +1,16 @@
+from spyro.meshing.meshing_functions import cells_per_wavelength
+import numpy as np
+
+
+def test_cpw_for_acoustic():
+    method = 'MLT'
+    degree = 3
+    dimension = 2
+    mlt3tri = cells_per_wavelength(method, degree, dimension)
+    dimension = 3
+    mlt3tet = cells_per_wavelength(method, degree, dimension)
+    assert np.isclose(mlt3tri, 3.70) and np.isclose(mlt3tet, 3.72)
+
+
+if __name__ == "__main__":
+    test_cpw_for_acoustic()

test/test_seismicmesh_integration.py

@@ -1,6 +1,15 @@
 import spyro
 import firedrake as fire
 import numpy as np
+import pytest
+
+
+def is_seismicmesh_installed():
+    try:
+        import SeismicMesh  # noqa: F401
+        return True
+    except ImportError:
+        return False
 
 
 def mean_edge_length(triangle):
@@ -14,6 +23,7 @@ def mean_edge_length(triangle):
     return (l0+l1+l2)/3.0
 
 
+@pytest.mark.skipif(not is_seismicmesh_installed(), reason="SeismicMesh is not installed")
 def test_spyro_seimicmesh_2d_homogeneous_generation():
     Lz = 1.0
     Lx = 2.0

test/test_velocity_smoother.py

@@ -0,0 +1,73 @@
+import firedrake as fire
+import spyro
+import segyio
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def get_vp_from_2dsegy(filename):
+    """
+    Extracts velocity profile (vp) data from a 2D SEG-Y file.
+    Parameters:
+    filename (str): The path to the SEG-Y file.
+    Returns:
+    np.ndarray: A 2D numpy array containing the velocity profile data,
+                with shape (nz, nx) where nz is the number of samples
+                and nx is the number of traces.
+    """
+
+    with segyio.open(filename, ignore_geometry=True) as f:
+        nz, nx = len(f.samples), len(f.trace)
+        show_vp = np.zeros(shape=(nz, nx))
+        for index, trace in enumerate(f.trace):
+            show_vp[:, index] = trace
+
+    return show_vp
+
+
+def test_write_segy_and_smooth(show=False):
+    vp_name = "velocity_models/test"
+    segy_file = vp_name + ".segy"
+    smoothed_file = "smoothed_test.segy"
+    mesh = fire.UnitSquareMesh(50, 50)
+    mesh.coordinates.dat.data[:, 0] *= -1
+
+    V = fire.FunctionSpace(mesh, "CG", 3)
+    x, y = fire.SpatialCoordinate(mesh)
+    r = 0.2
+    xc = -0.5
+    yc = 0.5
+
+    vp = fire.Function(V)
+
+    c = fire.conditional((x - xc) ** 2 + (y - yc) ** 2 < r**2, 3.0, 1.5)
+
+    vp.interpolate(c)
+
+    xi, yi, zi = spyro.io.write_function_to_grid(vp, V, 10.0 / 1000.0)
+    spyro.io.create_segy(zi, segy_file)
+    original_vp = get_vp_from_2dsegy(segy_file)
+
+    if show is True:
+        fig, ax = plt.subplots()
+        plt.pcolormesh(original_vp, shading="auto")
+        plt.title("Non smoothed model model")
+        plt.colorbar(label="P-wave velocity (km/s)")
+        plt.xlabel("x-direction (m)")
+        plt.ylabel("z-direction (m)")
+        ax.axis("equal")
+        plt.savefig("nonsmoothedtest.png")
+        plt.show()
+
+    spyro.tools.smooth_velocity_field_file(segy_file, smoothed_file, 5, show=show)
+
+    smoothed_vp = get_vp_from_2dsegy(smoothed_file)
+    check_boundary = np.isclose(original_vp[0, 0], smoothed_vp[0, 0])
+    check_centre = np.isclose(original_vp[48, 48], smoothed_vp[48, 48], rtol=1e-3)
+    check_halfway = original_vp[0, 0]*1.1 < smoothed_vp[24, 48] < original_vp[48, 48]*0.9
+
+    assert all([check_boundary, check_halfway, check_centre])
+
+
+if __name__ == "__main__":
+    test_write_segy_and_smooth(show=True)