Issue 0127 moment source

spyro/examples/camembert_elastic.py

@@ -24,7 +24,7 @@
 
 time_step = 2e-4  # [s]
 final_time = 1.5  # [s]
-out_freq = int(0.01/time_step)
+out_freq = int(0.1/time_step)
 
 n = 20
 mesh = fire.RectangleMesh(n, n, 0, L, originX=-L, diagonal='crossed')
@@ -57,7 +57,7 @@
     "source_type": "ricker",
     "source_locations": source_locations,
     "frequency": freq,
-    "delay": 0,
+    "delay": 1/freq,
     "delay_type": "time",
     "amplitude": np.array([0, smag]),
     # "amplitude": smag * np.eye(2),
@@ -78,7 +78,7 @@
     "final_time": final_time,
     "dt": time_step,
     "output_frequency": out_freq,
-    "gradient_sampling_frequency": 1,
+    "gradient_sampling_frequency": out_freq,
 }
 
 d["visualization"] = {

spyro/examples/elastic_analytical.py

@@ -0,0 +1,206 @@
+'''
+The analytical solutions are provided by the book:
+Quantitative Seismology (2nd edition) from Aki and Richards
+'''
+import numpy as np
+import os
+import spyro
+import time
+
+from firedrake.petsc import PETSc
+from math import pi as PI
+from mpi4py import MPI
+from numpy.linalg import norm
+from scipy.integrate import quad
+
+from spyro.utils.file_utils import mkdir_timestamp
+
+opts = PETSc.Options()
+
+moment_tensor = opts.getBool("moment_tensor", False)
+
+L = opts.getReal("L", default=450)  # Length (m)
+N = opts.getInt("N", default=30)    # Number of elements in each direction
+h = L/N                             # Element size (m)
+
+alpha = opts.getReal("alpha", default=1500)  # P-wave velocity
+beta = opts.getReal("beta", default=1000)    # S-wave velocity
+rho = opts.getReal("rho", default=2000)      # Density (kg/m3)
+
+smag = opts.getReal("amplitude", default=1e3)  # Source amplitude
+f0 = opts.getReal("frequency", default=20)     # Frequency (Hz)
+t0 = 1/f0                                      # Time delay
+
+tn = opts.getReal("total_time", default=0.3)  # Simulation time (s)
+nt = opts.getInt("time_steps", default=750)   # Number of time steps
+time_step = (tn/nt)
+out_freq = int(0.01/time_step)
+
+x_s = np.r_[-L/2, L/2, L/2]   # Source location (m)
+x_r = np.r_[opts.getReal("receiver_x", default=100),  # Receiver relative location (m)
+            opts.getReal("receiver_y", default=0),
+            opts.getReal("receiver_z", default=100)]
+
+
+def analytical_solution(i, j):
+    t = np.linspace(0, tn, nt)
+
+    u_near = np.zeros(nt)  # near field contribution
+    P_far = np.zeros(nt)   # P-wave far-field
+    S_far = np.zeros(nt)   # S-wave far field
+
+    r = np.linalg.norm(x_r)
+    gamma_i = x_r[i]/r
+    gamma_j = x_r[j]/r
+    delta_ij = 1 if i == j else 0
+
+    def X0(t):
+        a = PI*f0*(t - t0)
+        return (1 - 2*a**2)*np.exp(-a**2)
+
+    for k in range(nt):
+        res = quad(lambda tau: tau*X0(t[k] - tau), r/alpha, r/beta)
+        u_near[k] = smag*(1./(4*PI*rho))*(3*gamma_i*gamma_j - delta_ij)*(1./r**3)*res[0]
+        P_far[k] = smag*(1./(4*PI*rho*alpha**2))*gamma_i*gamma_j*(1./r)*X0(t[k] - r/alpha)
+        S_far[k] = smag*(1./(4*PI*rho*beta**2))*(gamma_i*gamma_j - delta_ij)*(1./r)*X0(t[k] - r/beta)
+
+    return u_near + P_far - S_far
+
+
+def explosive_source_analytical(i):
+    t = np.linspace(0, tn, nt)
+
+    P_mid = np.zeros(nt)  # P wave intermediate field
+    P_far = np.zeros(nt)  # P wave far field
+
+    r = np.linalg.norm(x_r)
+    gamma_i = x_r[i]/r
+
+    def w(t):
+        a = PI*f0*(t - t0)
+        return (t - t0)*np.exp(-a**2)
+
+    def w_dot(t):
+        a = PI*f0*(t - t0)
+        return (1 - 2*a**2)*np.exp(-a**2)
+
+    for k in range(nt):
+        P_mid[k] = smag*(gamma_i/(4*PI*rho*alpha**2))*(1./r**2)*w(t[k] - r/alpha)
+        P_far[k] = smag*(gamma_i/(4*PI*rho*alpha**3))*(1./r)*w_dot(t[k] - r/alpha)
+
+    return P_mid + P_far
+
+
+def numerical_solution(j):
+    if moment_tensor:
+        A0 = smag*np.eye(3)
+    else:
+        A0 = np.zeros(3)
+        A0[j] = smag
+
+    d = {}
+
+    d["options"] = {
+        "cell_type": "T",
+        "variant": "lumped",
+        "degree": 3,
+        "dimension": 3,
+    }
+
+    d["parallelism"] = {
+        "type": "automatic",
+    }
+
+    d["mesh"] = {
+        "Lz": L,
+        "Lx": L,
+        "Ly": L,
+        "mesh_file": None,
+        "mesh_type": "firedrake_mesh",
+    }
+
+    d["acquisition"] = {
+        "source_type": "ricker",
+        "source_locations": [x_s.tolist()],
+        "frequency": f0,
+        "delay": t0,
+        "delay_type": "time",
+        "amplitude": A0,
+        "receiver_locations": [(x_s + x_r).tolist()],
+    }
+
+    d["synthetic_data"] = {
+        "type": "object",
+        "density": rho,
+        "p_wave_velocity": alpha,
+        "s_wave_velocity": beta,
+        "real_velocity_file": None,
+    }
+
+    d["time_axis"] = {
+        "initial_time": 0,
+        "final_time": tn,
+        "dt": time_step,
+        "output_frequency": out_freq,
+        "gradient_sampling_frequency": 1,
+    }
+
+    d["visualization"] = {
+        "forward_output": True,
+        "forward_output_filename": "results/elastic_analytical.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={'edge_length': h})
+    wave.forward_solve()
+
+    return wave.receivers_output.reshape(nt + 1, 3)[0:-1, :]
+
+
+def err(u_a, u_n):
+    return norm(u_a - u_n)/norm(u_a)
+
+
+if __name__ == "__main__":
+    rank = MPI.COMM_WORLD.Get_rank()
+    if rank == 0:
+        start = time.time()
+
+    j = 0
+    U_x = analytical_solution(0, j)
+    U_y = analytical_solution(1, j)
+    U_z = analytical_solution(2, j)
+
+    u_n = numerical_solution(j)
+    u_x = u_n[:, 0]
+    u_y = u_n[:, 1]
+    u_z = u_n[:, 2]
+
+    if rank == 0:
+        end = time.time()
+        print(f"Elapsed time: {end - start:.2f} s")
+        print(f"err_x = {err(U_x, u_x)}")
+        print(f"err_y = {err(U_y, u_y)}")
+        print(f"err_z = {err(U_z, u_z)}")
+
+
+        if moment_tensor:
+            basename = "ExplosiveSource"
+        else:
+            basename = "ForceSource"
+        output_dir = mkdir_timestamp(basename)
+        np.save(os.path.join(output_dir, "x_analytical.npy"), U_x)
+        np.save(os.path.join(output_dir, "y_analytical.npy"), U_y)
+        np.save(os.path.join(output_dir, "z_analytical.npy"), U_z)
+        np.save(os.path.join(output_dir, "numerical.npy"), u_n)
+
+        print(f"Data saved to: {output_dir}")

spyro/io/basicio.py

@@ -53,9 +53,10 @@ def ensemble_plot(func):
     def wrapper(*args, **kwargs):
         num = args[0].number_of_sources
         _comm = args[0].comm
+        basename = kwargs.get("file_name", "")
         for snum in range(num):
             if is_owner(_comm, snum) and _comm.comm.rank == 0:
-                func(*args, **dict(kwargs, file_name=str(snum + 1)))
+                func(*args, **dict(kwargs, file_name=basename + str(snum + 1)))
 
     return wrapper
 

spyro/meshing/meshing_functions.py

@@ -356,8 +356,7 @@
         if self.dimension == 2:
             return self.create_seimicmesh_2d_mesh()
         elif self.dimension == 3:
-            raise NotImplementedError("Not implemented yet")
-            # return self.create_seismicmesh_3D_mesh()
+            return self.create_seismicmesh_3d_mesh()
         else:
             raise ValueError("dimension is not supported")
 
@@ -479,27 +478,67 @@
         )
 
         return fire.Mesh(self.output_file_name)
-        # raise NotImplementedError("Not implemented yet")
-
-
-# def create_firedrake_3D_mesh_based_on_parameters(dx, cell_type):
-#     nx = int(self.length_x / dx)
-#     nz = int(self.length_z / dx)
-#     ny = int(self.length_y / dx)
-#     if self.cell_type == "quadrilateral":
-#         quadrilateral = True
-#     else:
-#         quadrilateral = False
-
-#     return spyro.BoxMesh(
-#         nz,
-#         nx,
-#         ny,
-#         self.length_z,
-#         self.length_x,
-#         self.length_y,
-#         quadrilateral=quadrilateral,
-#     )
+
+    def create_seismicmesh_3d_mesh(self):
+        """
+        Creates a 3D mesh based on SeismicMesh meshing utilities.
+        """
+        if self.velocity_model is None:
+            return self.create_seismicmesh_3D_mesh_homogeneous()
+        else:
+            raise NotImplementedError("Not implemented yet")
+
+    def create_seismicmesh_3D_mesh_homogeneous(self):
+        """
+        Creates a 3D mesh based on SeismicMesh meshing utilities, with homogeneous velocity model.
+        """
+        Lz = self.length_z
+        Lx = self.length_x
+        Ly = self.length_y
+        pad = self.abc_pad
+
+        real_lz = Lz + pad
+        real_lx = Lx + 2 * pad
+        real_ly = Ly + 2 * pad
+
+        bbox = (-real_lz, 0.0,
+                -pad, real_lx - pad,
+                -pad, real_ly - pad)
+        cube = SeismicMesh.Cube(bbox)
+
+        points, cells = SeismicMesh.generate_mesh(
+            domain=cube,
+            edge_length=self.edge_length,
+            verbose=0,
+            max_iter=75,
+        )
+
+        points, cells = SeismicMesh.sliver_removal(
+            points=points,
+            domain=cube,
+            edge_length=self.edge_length,
+            preserve=True,
+        )
+
+        meshio.write_points_cells(
+            self.output_file_name,
+            points,
+            [("tetra", cells)],
+            file_format="gmsh22",
+            binary=False,
+        )
+
+        meshio.write_points_cells(
+            self.output_file_name + ".vtk",
+            points,
+            [("tetra", cells)],
+            file_format="vtk",
+        )
+
+        return fire.Mesh(self.output_file_name,
+                         distribution_parameters={
+                            "overlap_type": (fire.DistributedMeshOverlapType.NONE, 0)
+                        })
 
 
 def RectangleMesh(nx, ny, Lx, Ly, pad=None, comm=None, quadrilateral=False):

spyro/plots/plots.py

@@ -13,7 +13,7 @@
 def plot_shots(
     Wave_object,
     show=False,
-    file_name="1",
+    file_name="",
     vmin=-1e-5,
     vmax=1e-5,
     contour_lines=700,

spyro/receivers/dirac_delta_projector.py

@@ -492,10 +492,8 @@ def get_hexa_real_cell_node_map(V, mesh):
     cell_node_map = np.zeros(
         (layers * cells_per_layer, nodes_per_cell), dtype=int
     )
-    print(f"cnm size : {np.shape(cell_node_map)}", flush=True)
 
     for layer in range(layers):
-        print(f"layer : {layer} of {layers}", flush=True)
         for cell in range(cells_per_layer):
             cnm_base = weird_cnm[cell]
             cell_id = layer + layers * cell

spyro/solvers/acoustic_wave.py

@@ -210,3 +210,10 @@ def rhs_no_pml(self):
             return self.B.sub(0)
         else:
             return self.B
+
+    @override
+    def check_stability(self):
+        assert (
+            fire.norm(self.get_function()) < 1
+        ), "Numerical instability. Try reducing dt or building the " \
+           "mesh differently"

spyro/solvers/elastic_wave/isotropic_wave.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from firedrake import (assemble, Constant, curl, DirichletBC, div, Function,
+from firedrake import (assemble, norm, Constant, curl, DirichletBC, div, Function,
                        FunctionSpace, project, VectorFunctionSpace)
 
 from .elastic_wave import ElasticWave
@@ -234,3 +234,10 @@ def update_s_wave(self):
         self.s_wave.assign(project(curl(self.get_function()), self.C_h))
 
         return self.s_wave
+
+    @override
+    def check_stability(self):
+        assert (
+            np.isfinite(norm(self.get_function()))
+        ), "Numerical instability. Try reducing dt or building the " \
+           "mesh differently"