Merge pull request #6 from Olender/updating-readme

Updating readme
NDF-Poli-USP · Jul 25, 2023 · 9cd0aca · 9cd0aca
2 parents 77a13d8 + 16cf7fc
commit 9cd0aca
Show file tree

Hide file tree

Showing 6 changed files with 149 additions and 110 deletions.
diff --git a/README.md b/README.md
@@ -51,137 +51,119 @@ See the demos folder for an FWI example (this requires some other dependencies p
 ![Above shows the simulation at two timesteps in ParaView that results from running the code below](https://user-images.githubusercontent.com/18619644/94087976-7e81df00-fde5-11ea-96c0-474348286091.png)
 
 ```python
-from firedrake import (
-    RectangleMesh,
-    FunctionSpace,
-    Function,
-    SpatialCoordinate,
-    conditional,
-    File,
-)
-
 import spyro
 
-model = {}
-
-# Choose method and parameters
-model["opts"] = {
-    "method": "KMV",  # either CG or KMV
-    "quadratrue": "KMV", # Equi or KMV
-    "degree": 4,  # p order
-    "dimension": 2,  # dimension
+dictionary = {}
+
+# Choose spatial discretization method and parameters
+dictionary["options"] = {
+    # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q)
+    "cell_type": "T",  
+    # lumped, equispaced or DG, default is lumped "method":"MLT",
+    # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral)
+    # You can either specify a cell_type+variant or a method.
+    "variant": 'lumped',  
+    # Polynomial order of the spatial discretion's basis functions.
+    # For MLT we recomend 4th order in 2D, 3rd order in 3D, for SEM 4th or 8th.
+    "degree": 4,  
+    # Dimension (2 or 3)
+    "dimension": 2,  
 }
 
-# Number of cores for the shot. For simplicity, we keep things serial.
-# spyro however supports both spatial parallelism and "shot" parallelism.
-model["parallelism"] = {
-    "type": "automatic",  # options: automatic (same number of cores for evey processor) or spatial
+# Number of cores for the shot. For simplicity, we keep things automatic.
+# spyro supports both spatial parallelism and "shot" parallelism.
+dictionary["parallelism"] = {
+    # options: automatic (same number of cores for every shot) or spatial
+    "type": "automatic",
 }
 
 # Define the domain size without the PML. Here we'll assume a 0.75 x 1.50 km
-# domain and reserve the remaining 250 m for the Perfectly Matched Layer (PML) to absorb
-# outgoing waves on three sides (eg., -z, +-x sides) of the domain.
-model["mesh"] = {
-    "Lz": 0.75,  # depth in km - always positive
-    "Lx": 1.5,  # width in km - always positive
-    "Ly": 0.0,  # thickness in km - always positive
-    "meshfile": "not_used.msh",
-    "initmodel": "not_used.hdf5",
-    "truemodel": "not_used.hdf5",
-}
-
-# Specify a 250-m PML on the three sides of the domain to damp outgoing waves.
-model["BCs"] = {
-    "status": True,  # True or false
-    "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.7,  # maximum acoustic wave velocity in PML - km/s
-    "R": 1e-6,  # theoretical reflection coefficient
-    "lz": 0.25,  # thickness of the PML in the z-direction (km) - always positive
-    "lx": 0.25,  # 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
+dictionary["mesh"] = {
+    # depth in km - always positive
+    "Lz": 0.75,
+    # width in km - always positive
+    "Lx": 1.50,
+    # thickness in km - always positive
+    "Ly": 0.0,
+    # If we are loading and external .msh mesh file
+    "mesh_file": None,
+    # options: None (default), firedrake_mesh, user_mesh, or SeismicMesh
+    # use this opion if your are not loading an external file
+    # 'firedrake_mesh' will create an automatic mesh using firedrake's built-in meshing tools
+    # 'user_mesh' gives the option to load other user generated meshes from unsuported formats
+    # 'SeismicMesh' automatically creates a waveform adapted unstructured mesh to reduce total
+    # DoFs using the SeismicMesh tool.
+    "mesh_type": "firedrake_mesh",
 }
 
 # Create a source injection operator. Here we use a single source with a
 # Ricker wavelet that has a peak frequency of 8 Hz injected at the center of the mesh.
 # We also specify to record the solution at 101 microphones near the top of the domain.
 # This transect of receivers is created with the helper function `create_transect`.
-model["acquisition"] = {
-    "source_type": "Ricker",
-    "source_pos": [(-0.1, 0.75)],
+dictionary["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": [(-0.3, 0.75)],
     "frequency": 8.0,
     "delay": 1.0,
     "receiver_locations": spyro.create_transect(
-        (-0.10, 0.1), (-0.10, 1.4), 100
+        (-0.5, 0.1), (-0.5, 1.4), 100
     ),
 }
 
 # Simulate for 2.0 seconds.
-model["timeaxis"] = {
-    "t0": 0.0,  #  Initial time for event
-    "tf": 2.00,  # Final time for event
-    "dt": 0.0005,  # timestep size
-    "amplitude": 1,  # the Ricker has an amplitude of 1.
-    "nspool": 100,  # how frequently to output solution to pvds
-    "fspool": 100,  # how frequently to save solution to RAM
+dictionary["time_axis"] = {
+    # Initial time for event
+    "initial_time": 0.0,
+    # Final time for event
+    "final_time": 0.50,
+    # timestep size
+    "dt": 0.0001,
+    # the Ricker has an amplitude of 1.
+    "amplitude": 1,
+    # how frequently to output solution to pvds
+    "output_frequency": 100,
+    # how frequently to save solution to RAM
+    "gradient_sampling_frequency": 100,
 }
 
+dictionary["visualization"] = {
+    "forward_output" : True,
+    "output_filename": "results/forward_output.pvd",
+    "fwi_velocity_model_output": False,
+    "velocity_model_filename": None,
+    "gradient_output": False,
+    "gradient_filename": None,
+}
 
-# Create a simple mesh of a rectangle ∈ [1 x 2] km with ~100 m sized elements
-# and then create a function space for P=1 Continuous Galerkin FEM
-mesh = RectangleMesh(100, 200, 1.0, 2.0)
-
-# We edit the coordinates of the mesh so that it's in the (z, x) plane
-# and has a domain padding of 250 m on three sides, which will be used later to show
-# the Perfectly Matched Layer (PML). More complex 2D/3D meshes can be automatically generated with
-# SeismicMesh https://github.com/krober10nd/SeismicMesh
-mesh.coordinates.dat.data[:, 0] -= 1.0
-mesh.coordinates.dat.data[:, 1] -= 0.25
 
+# Create an AcousticWave object with the above dictionary.
+Wave_obj = spyro.AcousticWave(dictionary=dictionary)
 
-# Create the computational environment
-comm = spyro.utils.mpi_init(model)
+# Defines the element size in the automatically generated firedrake mesh.
+Wave_obj.set_mesh(dx=0.01)
 
-element = spyro.domains.space.FE_method(
-    mesh, model["opts"]["method"], model["opts"]["degree"]
-)
-V = FunctionSpace(mesh, element)
 
-# Manually create a simple two layer seismic velocity model `vp`.
-# Note: the user can specify their own velocity model in a HDF5 file format
-# in the above two lines using SeismicMesh.
-# If so, the HDF5 file has to contain an array with
+# Manually create a simple two layer seismic velocity model.
+# Note: the user can specify their own velocity model in a HDF5 or SEG-Y file format.
+# The HDF5 file has to contain an array with
 # the velocity data and it is linearly interpolated onto the mesh nodes at run-time.
-x, y = SpatialCoordinate(mesh)
-velocity = conditional(x > -0.35, 1.5, 3.0)
-vp = Function(V, name="velocity").interpolate(velocity)
-# These pvd files can be easily visualized in ParaView!
-File("simple_velocity_model.pvd").write(vp)
-
-
-# Now we instantiate both the receivers and source objects.
-sources = spyro.Sources(model, mesh, V, comm)
-
-receivers = spyro.Receivers(model, mesh, V, comm)
-
-# Create a wavelet to force the simulation
-wavelet = spyro.full_ricker_wavelet(dt=0.0005, tf=2.0, freq=8.0)
+z = Wave_obj.mesh_z
+import firedrake as fire
+velocity_conditional = fire.conditional(z > -0.35, 1.5, 3.0)
+Wave_obj.set_initial_velocity_model(conditional=velocity_conditional, output=True)
 
 # And now we simulate the shot using a 2nd order central time-stepping scheme
 # Note: simulation results are stored in the folder `~/results/` by default
-p_field, p_at_recv = spyro.solvers.forward(
-    model, mesh, comm, vp, sources, wavelet, receivers
-)
+Wave_obj.forward_solve()
 
 # Visualize the shot record
-spyro.plots.plot_shots(model, comm, p_at_recv)
+spyro.plots.plot_shots(Wave_obj, show=True)
 
 # Save the shot (a Numpy array) as a pickle for other use.
-spyro.io.save_shots(model, comm, p_at_recv)
+spyro.io.save_shots(Wave_obj)
 
 # can be loaded back via
-my_shot = spyro.io.load_shots(model, comm)
+my_shot = spyro.io.load_shots(Wave_obj)
 ```
 
 ### Testing

diff --git a/spyro/io/basicio.py b/spyro/io/basicio.py
@@ -50,9 +50,8 @@ def ensemble_plot(func):
     ensemble parallelism"""
 
     def wrapper(*args, **kwargs):
-        acq = args[0].get("acquisition")
-        num = len(acq["source_pos"])
-        _comm = args[1]
+        num = args[0].number_of_sources
+        _comm = args[0].comm
         for snum in range(num):
             if is_owner(_comm, snum) and _comm.comm.rank == 0:
                 func(*args, **dict(kwargs, file_name=str(snum + 1)))

diff --git a/spyro/plots/plots.py b/spyro/plots/plots.py
@@ -8,9 +8,7 @@
 
 @ensemble_plot
 def plot_shots(
-    model,
-    comm,
-    arr,
+    Wave_object,
     show=False,
     file_name="1",
     vmin=-1e-5,
@@ -49,12 +47,14 @@ def plot_shots(
     None
     """
 
-    num_recvs = len(model["acquisition"]["receiver_locations"])
+    num_recvs = Wave_object.number_of_receivers
 
-    dt = model["timeaxis"]["dt"]
-    tf = model["timeaxis"]["tf"]
+    dt = Wave_object.dt
+    tf = Wave_object.final_time
 
-    nt = int(tf / dt)  # number of timesteps
+    arr = Wave_object.receivers_output
+
+    nt = int(tf / dt)  + 1 # number of timesteps
 
     if end_index == 0:
         end_index = num_recvs

diff --git a/spyro/solvers/wave.py b/spyro/solvers/wave.py
@@ -118,6 +118,7 @@ def set_initial_velocity_model(
         velocity_model_function=None,
         expression=None,
         new_file=None,
+        output=False,
     ):
         """Method to define new user velocity model or file. It is optional.
 
@@ -151,9 +152,6 @@ def set_initial_velocity_model(
             V = self.function_space
             vp = fire.Function(V, name="velocity")
             vp.interpolate(expression)
-            fire.File("initial_velocity_model_3d.pvd").write(
-                vp, name="velocity"
-            )
             self.initial_velocity_model = vp
         elif velocity_model_function is not None:
             self.initial_velocity_model = velocity_model_function
@@ -163,13 +161,16 @@ def set_initial_velocity_model(
             V = self.function_space
             vp = fire.Function(V, name="velocity")
             vp.interpolate(fire.Constant(constant))
-            fire.File("initial_velocity_model.pvd").write(vp, name="velocity")
             self.initial_velocity_model = vp
         else:
             raise ValueError(
                 "Please specify either a conditional, expression, firedrake \
                     function or new file name (segy or hdf5)."
             )
+        if output:
+            fire.File("initial_velocity_model.pvd").write(
+                self.initial_velocity_model, name="velocity"
+            )
 
     def _get_initial_velocity_model(self):
         if self.velocity_model_type == "conditional":

diff --git a/test/inputfiles/model.py b/test/inputfiles/model.py
@@ -0,0 +1,57 @@
+
+
+dictionary = {}
+dictionary["options"] = {
+    "cell_type": "Q",  # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q)
+    "variant": 'lumped',  # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method
+    "degree": 4,  # p order
+    "dimension": 2,  # dimension
+}
+
+# Number of cores for the shot. For simplicity, we keep things serial.
+# spyro however supports both spatial parallelism and "shot" parallelism.
+dictionary["parallelism"] = {
+    "type": "automatic",  # options: automatic (same number of cores for evey processor) or spatial
+}
+
+# Define the domain size without the PML. Here we'll assume a 1.00 x 1.00 km
+# domain and reserve the remaining 250 m for the Perfectly Matched Layer (PML) to absorb
+# outgoing waves on three sides (eg., -z, +-x sides) of the domain.
+dictionary["mesh"] = {
+    "Lz": 1.0,  # depth in km - always positive
+    "Lx": 1.0,  # width in km - always positive
+    "Ly": 0.0,  # thickness in km - always positive
+    "mesh_type": "firedrake_mesh",  # options: firedrake_mesh or user_mesh
+    "mesh_file": None,  # specify the mesh file
+}
+
+# Create a source injection operator. Here we use a single source with a
+# Ricker wavelet that has a peak frequency of 5 Hz injected at the center of the mesh.
+# We also specify to record the solution at a microphone near the top of the domain.
+# This transect of receivers is created with the helper function `create_transect`.
+dictionary["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": [(-1.0, 1.0)],#, (-0.605, 1.7), (-0.61, 1.7), (-0.615, 1.7)],#, (-0.1, 1.5), (-0.1, 2.0), (-0.1, 2.5), (-0.1, 3.0)],
+    "frequency": 5.0,
+    "delay": 1.5,
+    "receiver_locations": [(-0.0, 0.5)],
+}
+
+# Simulate for 2.0 seconds.
+dictionary["time_axis"] = {
+    "initial_time": 0.0,  # Initial time for event
+    "final_time": 1.0,  # Final time for event
+    "dt": 0.0005,  # timestep size
+    "amplitude": 1,  # the Ricker has an amplitude of 1.
+    "output_frequency": 100,  # how frequently to output solution to pvds
+    "gradient_sampling_frequency": 1,  # how frequently to save solution to RAM
+}
+
+dictionary["visualization"] = {
+    "forward_output" : False,
+    "output_filename": "results/forward_output.pvd",
+    "fwi_velocity_model_output": False,
+    "velocity_model_filename": None,
+    "gradient_output": False,
+    "gradient_filename": None,
+}
diff --git a/test/test_io.py b/test/test_io.py
@@ -83,7 +83,7 @@ def test_read_and_write_segy():
 
 
 def test_saving_shot_record():
-    from model import dictionary
+    from .inputfiles.model import dictionary
     dictionary["time_axis"]["final_time"] = 0.5
     Wave_obj = spyro.AcousticWave(dictionary=dictionary)
     Wave_obj.set_mesh(dx=0.02)
@@ -93,7 +93,7 @@ def test_saving_shot_record():
 
 
 def test_loading_shot_record():
-    from model import dictionary
+    from .inputfiles.model import dictionary
     dictionary["time_axis"]["final_time"] = 0.5
     Wave_obj = spyro.AcousticWave(dictionary=dictionary)
     Wave_obj.set_mesh(dx=0.02)