Merge pull request #84 from NDF-Poli-USP/switching-to-soft-dependency…

…-for-seismicmesh Switching to soft dependency for seismicmesh
NDF-Poli-USP · Jan 17, 2023 · 6f14329 · 6f14329
2 parents 27f87b1 + 139ac2a
commit 6f14329
Show file tree

Hide file tree

Showing 60 changed files with 1,465 additions and 1,135 deletions.
diff --git a/README.md b/README.md
@@ -14,6 +14,8 @@ To use Spyro, you'll need to have some knowledge of Python and some basic concep
 
 Discussions about development take place on our Slack channel. Everyone is invited to join using the link: https://join.slack.com/t/spyroworkspace/shared_invite/zt-u87ih28m-2h9JobfkdArs4ku3a1wLLQ
 
+If you want to know more or cite our code please see our open access publication: https://gmd.copernicus.org/articles/15/8639/2022/gmd-15-8639-2022.html
+
 Functionality
 =============
 

diff --git a/benchmarks/benchmark_2d.py b/benchmarks/benchmark_2d.py
@@ -1,10 +1,13 @@
-import SeismicMesh
 import meshio
 import firedrake as fire
 import spyro
 import time
+try:
+    import SeismicMesh
+except ImportError:
+    pass
 
-## Adding model parameters:
+# Adding model parameters:
 model = {}
 
 model["opts"] = {
@@ -16,7 +19,8 @@
 
 model["parallelism"] = {
     "type": "spatial",  # options: automatic, spatial, or custom.
-    "custom_cores_per_shot": [],  # only if the user wants a different number of cores for every shot.
+    "custom_cores_per_shot": [],  # only if the user wants a different number
+    # of cores for every shot.
     # input is a list of integers with the length of the number of shots.
 }
 
@@ -31,14 +35,15 @@
 
 model["BCs"] = {
     "status": True,  # True or false
-    "outer_bc": "non-reflective",  #  None or non-reflective (outer boundary condition)
+    # None or non-reflective (outer boundary condition)
+    "outer_bc": "non-reflective",
     "damping_type": "polynomial",  # polynomial, hyperbolic, shifted_hyperbolic
     "exponent": 2,  # damping layer has a exponent variation
     "cmax": 4.5,  # maximum acoustic wave velocity in PML - km/s
     "R": 1e-6,  # theoretical reflection coefficient
-    "lz": 0.286,  # thickness of the PML in the z-direction (km) - always positive
-    "lx": 0.286,  # thickness of the PML in the x-direction (km) - always positive
-    "ly": 0.0,  # thickness of the PML in the y-direction (km) - always positive
+    "lz": 0.286,  # thickness of the PML in z-direction (km)-always positive
+    "lx": 0.286,  # thickness of the PML in x-direction (km)-always positive
+    "ly": 0.0,  # thickness of the PML in y-direction (km)-always positive
 }
 
 model["acquisition"] = {
@@ -49,14 +54,14 @@
     "delay": 1.0,
     "num_receivers": 15,
     "receiver_locations": spyro.create_transect(
-        (-5.72,  6.29),
+        (-5.72, 6.29),
         (-5.72, 14.29),
         15),
 }
 
 # Simulate for 1.0 seconds.
 model["timeaxis"] = {
-    "t0": 0.0,  #  Initial time for event
+    "t0": 0.0,  # Initial time for event
     "tf": 4.0,  # Final time for event
     "dt": 0.0005,  # timestep size
     "amplitude": 1,  # the Ricker has an amplitude of 1.
@@ -66,11 +71,12 @@
 
 comm = spyro.utils.mpi_init(model)
 
-## Starting meshing procedure with seismic mesh. This can be done seperately in order to not have to
+# Starting meshing procedure with seismic mesh. This can be done seperately
+# in order to not have to
 # when testing multiple cores.
 
-print('Entering mesh generation', flush = True)
-if model['opts']['degree']   == 2:
+print('Entering mesh generation', flush=True)
+if model['opts']['degree'] == 2:
     M = 5.85
 elif model['opts']['degree'] == 3:
     M = 3.08
@@ -84,25 +90,30 @@
 Lx = model["mesh"]["Lx"]
 pad = model["BCs"]["lz"]
 
-bbox = (-Real_Lz, 0.0,-pad, Lx+pad)
+bbox = (-Real_Lz, 0.0, -pad, Lx+pad)
 rec = SeismicMesh.Rectangle(bbox)
 if comm.comm.rank == 0:
     points, cells = SeismicMesh.generate_mesh(
-        domain=rec, 
-        edge_length=edge_length, 
-        comm = comm.ensemble_comm,
-        verbose = 0
-        )
-
-    points, cells = SeismicMesh.geometry.delete_boundary_entities(points, cells, min_qual= 0.6)
-
-    meshio.write_points_cells("meshes/benchmark_2d.msh",
-        points,[("triangle", cells)],
-        file_format="gmsh22", 
-        binary = False
-        )
-
-## Mesh generation finishes here.
+        domain=rec,
+        edge_length=edge_length,
+        comm=comm.ensemble_comm,
+        verbose=0
+    )
+
+    points, cells = SeismicMesh.geometry.delete_boundary_entities(
+        points,
+        cells,
+        min_qual=0.6,
+    )
+
+    meshio.write_points_cells(
+        "meshes/benchmark_2d.msh",
+        points, [("triangle", cells)],
+        file_format="gmsh22",
+        binary=False
+    )
+
+# Mesh generation finishes here.
 
 mesh = fire.Mesh(
     "meshes/benchmark_2d.msh",
@@ -117,7 +128,12 @@
 if comm.ensemble_comm.rank == 0 and comm.comm.rank == 0:
     print(f"Setting up {method} a {degree}tetra element", flush=True)
 
-element = fire.FiniteElement(method, mesh.ufl_cell(), degree=degree, variant="KMV")
+element = fire.FiniteElement(
+    method,
+    mesh.ufl_cell(),
+    degree=degree,
+    variant="KMV"
+)
 
 V = fire.FunctionSpace(mesh, element)
 
@@ -139,5 +155,13 @@
 )
 
 t1 = time.time()
-p, p_r = spyro.solvers.forward(model, mesh, comm, vp, sources, wavelet, receivers)
+p, p_r = spyro.solvers.forward(
+    model,
+    mesh,
+    comm,
+    vp,
+    sources,
+    wavelet,
+    receivers
+)
 print(time.time() - t1, flush=True)
diff --git a/benchmarks/benchmark_3d.py b/benchmarks/benchmark_3d.py
@@ -1,10 +1,13 @@
-import SeismicMesh
 import meshio
 import firedrake as fire
 import spyro
 import time
+try:
+    import SeismicMesh
+except ImportError:
+    pass
 
-## Adding model parameters:
+# Adding model parameters:
 model = {}
 
 model["opts"] = {
@@ -16,7 +19,8 @@
 
 model["parallelism"] = {
     "type": "spatial",  # options: automatic, spatial, or custom.
-    "custom_cores_per_shot": [],  # only if the user wants a different number of cores for every shot.
+    # only if the user wants a different number of cores for every shot.
+    "custom_cores_per_shot": [],
     # input is a list of integers with the length of the number of shots.
 }
 
@@ -31,14 +35,15 @@
 
 model["BCs"] = {
     "status": True,  # True or false
-    "outer_bc": "non-reflective",  #  None or non-reflective (outer boundary condition)
+    # None or non-reflective (outer boundary condition)
+    "outer_bc": "non-reflective",
     "damping_type": "polynomial",  # polynomial, hyperbolic, shifted_hyperbolic
     "exponent": 1,  # damping layer has a exponent variation
     "cmax": 4.7,  # maximum acoustic wave velocity in PML - km/s
     "R": 1e-3,  # theoretical reflection coefficient
-    "lz": 0.286,  # thickness of the PML in the z-direction (km) - always positive
-    "lx": 0.286,  # thickness of the PML in the x-direction (km) - always positive
-    "ly": 0.286,  # thickness of the PML in the y-direction (km) - always positive
+    "lz": 0.286,  # thickness of the PML in z-direction (km) - always positive
+    "lx": 0.286,  # thickness of the PML in x-direction (km) - always positive
+    "ly": 0.286,  # thickness of the PML in y-direction (km) - always positive
 }
 
 model["acquisition"] = {
@@ -56,7 +61,7 @@
 
 # Simulate for 1.0 seconds.
 model["timeaxis"] = {
-    "t0": 0.0,  #  Initial time for event
+    "t0": 0.0,  # Initial time for event
     "tf": 1.0,  # Final time for event
     "dt": 0.005,  # timestep size
     "amplitude": 1,  # the Ricker has an amplitude of 1.
@@ -66,11 +71,12 @@
 
 comm = spyro.utils.mpi_init(model)
 
-## Starting meshing procedure with seismic mesh. This can be done seperately in order to not have to
+# Starting meshing procedure with seismic mesh. This can be done seperately
+# in order to not have to
 # when testing multiple cores.
 
-print('Entering mesh generation', flush = True)
-if model['opts']['degree']   == 2:
+print('Entering mesh generation', flush=True)
+if model['opts']['degree'] == 2:
     M = 5.1
 elif model['opts']['degree'] == 3:
     M = 3.1
@@ -81,25 +87,34 @@
 Ly = model["mesh"]["Ly"]
 pad = model["BCs"]["lz"]
 
-bbox = (-Real_Lz, 0.0,-pad, Lx+pad, -pad, Ly+pad)
+bbox = (-Real_Lz, 0.0, -pad, Lx+pad, -pad, Ly+pad)
 cube = SeismicMesh.Cube(bbox)
 points, cells = SeismicMesh.generate_mesh(
-    domain=cube, 
-    edge_length=edge_length, 
-    max_iter = 80,
-    comm = comm.ensemble_comm,
-    verbose = 2
-    )
+    domain=cube,
+    edge_length=edge_length,
+    max_iter=80,
+    comm=comm.ensemble_comm,
+    verbose=2
+)
 
-points, cells = SeismicMesh.sliver_removal(points=points, bbox=bbox, max_iter=100, domain=cube, edge_length=edge_length, preserve=True)
+points, cells = SeismicMesh.sliver_removal(
+    points=points,
+    bbox=bbox,
+    max_iter=100,
+    domain=cube,
+    edge_length=edge_length,
+    preserve=True
+)
 
-meshio.write_points_cells("meshes/benchmark_3d.msh",
-    points,[("tetra", cells)],
-    file_format="gmsh22", 
-    binary = False
-    )
+meshio.write_points_cells(
+    "meshes/benchmark_3d.msh",
+    points,
+    [("tetra", cells)],
+    file_format="gmsh22",
+    binary=False
+)
 
-## Mesh generation finishes here.
+# Mesh generation finishes here.
 
 mesh = fire.Mesh(
     "meshes/benchmark_3d.msh",
@@ -114,7 +129,12 @@
 if comm.ensemble_comm.rank == 0 and comm.comm.rank == 0:
     print(f"Setting up {method} a {degree}tetra element", flush=True)
 
-element = fire.FiniteElement(method, mesh.ufl_cell(), degree=degree, variant="KMV")
+element = fire.FiniteElement(
+    method,
+    mesh.ufl_cell(),
+    degree=degree,
+    variant="KMV"
+)
 
 V = fire.FunctionSpace(mesh, element)
 
@@ -136,5 +156,13 @@
 )
 
 t1 = time.time()
-p, p_r = spyro.solvers.forward(model, mesh, comm, vp, sources, wavelet, receivers)
+p, p_r = spyro.solvers.forward(
+    model,
+    mesh,
+    comm,
+    vp,
+    sources,
+    wavelet,
+    receivers
+)
 print(time.time() - t1, flush=True)
diff --git a/benchmarks/benchmark_cuboid.py b/benchmarks/benchmark_cuboid.py
@@ -32,7 +32,7 @@
 
 model["PML"] = {
     "status": True,  # True or false
-    "outer_bc": "non-reflective",  #  None or non-reflective (outer boundary condition)
+    "outer_bc": "non-reflective",  # None or non-reflective (outer boundary condition)
     "damping_type": "polynomial",  # polynomial, hyperbolic, shifted_hyperbolic
     "exponent": 2,  # damping layer has a exponent variation
     "cmax": 1.0,  # maximum acoustic wave velocity in PML - km/s
@@ -42,6 +42,7 @@
     "ly": 1.0,  # thickness of the PML in the y-direction (km) - always positive
 }
 
+
 model["acquisition"] = {
     "source_type": "Ricker",
     "num_sources": 1,
@@ -56,7 +57,7 @@
 
 # Simulate for 1.0 seconds.
 model["timeaxis"] = {
-    "t0": 0.0,  #  Initial time for event
+    "t0": 0.0,  # Initial time for event
     "tf": 1.0,  # Final time for event
     "dt": 0.0005,  # timestep size
     "amplitude": 1,  # the Ricker has an amplitude of 1.

diff --git a/demos/run_forward.py b/demos/run_forward.py
@@ -21,7 +21,7 @@
 }
 model["BCs"] = {
     "status": True,  # True or false
-    "outer_bc": "non-reflective",  #  None or non-reflective (outer boundary condition)
+    "outer_bc": "non-reflective",  # None or non-reflective (outer boundary condition)
     "damping_type": "polynomial",  # polynomial, hyperbolic, shifted_hyperbolic
     "exponent": 2,  # damping layer has a exponent variation
     "cmax": 4.5,  # maximum acoustic wave velocity in PML - km/s
@@ -40,7 +40,7 @@
     "receiver_locations": spyro.create_transect((-0.10, 0.1), (-0.10, 17.0), 500),
 }
 model["timeaxis"] = {
-    "t0": 0.0,  #  Initial time for event
+    "t0": 0.0,  # Initial time for event
     "tf": 6.00,  # Final time for event
     "dt": 0.001,
     "amplitude": 1,  # the Ricker has an amplitude of 1.

diff --git a/demos/run_fwi.py b/demos/run_fwi.py
@@ -20,7 +20,7 @@
     "degree": 5,  # p order
     "dimension": 2,  # dimension
     "regularization": True,  # regularization is on?
-    "gamma": 1.e-6, # regularization parameter
+    "gamma": 1.0e-6,  # regularization parameter
 }
 model["parallelism"] = {
     "type": "automatic",
@@ -35,7 +35,7 @@
 }
 model["BCs"] = {
     "status": True,  # True or false
-    "outer_bc": "non-reflective",  #  None or non-reflective (outer boundary condition)
+    "outer_bc": "non-reflective",  # None or non-reflective (outer boundary condition)
     "damping_type": "polynomial",  # polynomial, hyperbolic, shifted_hyperbolic
     "exponent": 2,  # damping layer has a exponent variation
     "cmax": 4.5,  # maximum acoustic wave velocity in PML - km/s
@@ -54,7 +54,7 @@
     "receiver_locations": spyro.create_transect((-0.10, 0.1), (-0.10, 17.0), 500),
 }
 model["timeaxis"] = {
-    "t0": 0.0,  #  Initial time for event
+    "t0": 0.0,  # Initial time for event
     "tf": 6.00,  # Final time for event
     "dt": 0.001,
     "amplitude": 1,  # the Ricker has an amplitude of 1.
@@ -179,8 +179,8 @@ def gradient(self, g, x, tol):
         if comm.comm.size > 1:
             dJ /= comm.comm.size
         # regularize the gradient if asked.
-        if model['opts']['regularization']:
-            gamma = model['opts']['gamma']
+        if model["opts"]["regularization"]:
+            gamma = model["opts"]["gamma"]
             dJ = regularize_gradient(vp, dJ, gamma)
         # mask the water layer
         dJ.dat.data[water] = 0.0