switching shots

spyro/examples/__init__.py

@@ -1,4 +1,5 @@
 from .camembert import Camembert_acoustic
+from .camembert import Camembert_acoustic_FWI
 from .marmousi import Marmousi_acoustic
 from .cut_marmousi import Cut_marmousi_acoustic
 from .example_model import Example_model_acoustic
@@ -10,6 +11,7 @@
 
 __all__ = [
     "Camembert_acoustic",
+    "Camembert_acoustic_FWI",
     "Marmousi_acoustic",
     "Example_model_acoustic",
     "Example_model_acoustic_FWI",

spyro/examples/camembert.py

@@ -1,5 +1,5 @@
 from spyro import create_transect
-from spyro.examples.rectangle import Rectangle_acoustic
+from spyro.examples.rectangle import Rectangle_acoustic, Rectangle_acoustic_FWI
 import firedrake as fire
 
 camembert_optimization_parameters = {
@@ -135,11 +135,13 @@ def __init__(
         dictionary=None,
         example_dictionary=camembert_dictionary,
         comm=None,
+        periodic=False,
     ):
         super().__init__(
             dictionary=dictionary,
             example_dictionary=example_dictionary,
             comm=comm,
+            periodic=False,
         )
         self._camembert_velocity_model()
 
@@ -156,3 +158,57 @@ def _camembert_velocity_model(self):
         )
         self.set_initial_velocity_model(conditional=cond, dg_velocity_model=False)
         return None
+
+
+camembert_dictionary["inversion"] = {
+    "perform_fwi": True,  # switch to true to make a FWI
+    "initial_guess_model_file": None,
+    "shot_record_file": None,
+}
+
+
+class Camembert_acoustic_FWI(Rectangle_acoustic_FWI):
+    """Camembert model.
+    This class is a child of the Example_model class.
+    It is used to create a dictionary with the parameters of the
+    Camembert model.
+
+    Parameters
+    ----------
+    dictionary : dict, optional
+        Dictionary with the parameters of the model that are different from
+        the default Camembert model. The default is None.
+
+    """
+
+    def __init__(
+        self,
+        dictionary=None,
+        example_dictionary=camembert_dictionary,
+        comm=None,
+        periodic=False,
+    ):
+        super().__init__(
+            dictionary=dictionary,
+            example_dictionary=example_dictionary,
+            comm=comm,
+            periodic=False,
+        )
+        self._camembert_velocity_model()
+        self.real_velocity_model = self.initial_velocity_model
+        self.real_mesh = self.mesh
+
+    def _camembert_velocity_model(self):
+        camembert_dict = self.input_dictionary["camembert_options"]
+        z = self.mesh_z
+        x = self.mesh_x
+        zc, xc = camembert_dict["circle_center"]
+        rc = camembert_dict["radius"]
+        c_salt = camembert_dict["inside_circle_velocity"]
+        c_not_salt = camembert_dict["outside_velocity"]
+        cond = fire.conditional(
+            (z - zc) ** 2 + (x - xc) ** 2 < rc**2, c_salt, c_not_salt
+        )
+        self.set_initial_velocity_model(conditional=cond, dg_velocity_model=False)
+        return None
+

spyro/examples/immersed_polygon.py

@@ -40,8 +40,8 @@
     # every processor) or spatial
 }
 polygon_dictionary["mesh"] = {
-    "Lz": 2.0,  # depth in km - always positive
-    "Lx": 3.0,  # width in km - always positive
+    "Lz": 1.0,  # depth in km - always positive
+    "Lx": 1.0,  # width in km - always positive
     "Ly": 0.0,  # thickness in km - always positive
     "h": 0.05,  # mesh size in km
     "mesh_file": None,
@@ -149,7 +149,7 @@ def _polygon_velocity_model(self):
             cond = fire.conditional(z <= d1, v2, v1)
         cond = fire.conditional(z <= d2 - 0.2*x, vl, cond)
 
-        cond = fire.conditional(300*((x-1.5)*(-z-0.7))**2 + ((x-1.5)+(-z-0.7))**2 <= 0.300**2, v2+dv, cond)
+        cond = fire.conditional(300*((x-0.5)*(-z-0.5))**2 + ((x-0.5)+(-z-0.5))**2 <= 0.300**2, v2+dv, cond)
 
         if self.abc_pad_length > 0.0:
             middle_of_pad = -self.length_z - self.abc_pad_length*0.5

spyro/solvers/inversion.py

@@ -380,6 +380,7 @@ def get_functional(self, c=None):
 
         self.functional_history.append(Jm)
         self.functional = Jm
+        print(f"Functional: {Jm} at iteration: {self.current_iteration}", flush=True)
 
         return Jm
 

test_camembert_fwi.py

@@ -0,0 +1,100 @@
+# from mpi4py.MPI import COMM_WORLD
+# import debugpy
+# debugpy.listen(3000 + COMM_WORLD.rank)
+# debugpy.wait_for_client()
+import spyro
+import firedrake as fire
+import numpy as np
+import matplotlib.pyplot as plt
+import warnings
+warnings.filterwarnings("ignore")
+
+
+final_time = 0.9
+
+dictionary = {}
+dictionary["options"] = {
+    "cell_type": "T",  # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q)
+    "variant": "lumped",  # lumped, equispaced or DG, default is lumped
+    "degree": 4,  # p order
+    "dimension": 2,  # dimension
+}
+dictionary["parallelism"] = {
+    "type": "automatic",  # options: automatic (same number of cores for evey processor) or spatial
+}
+dictionary["mesh"] = {
+    "Lz": 2.0,  # depth in km - always positive   # Como ver isso sem ler a malha?
+    "Lx": 2.0,  # width in km - always positive
+    "Ly": 0.0,  # thickness in km - always positive
+    "mesh_file": None,
+    "mesh_type": "firedrake_mesh",
+}
+dictionary["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": spyro.create_transect((-0.55, 0.7), (-0.55, 1.3), 6),
+    # "source_locations": [(-1.1, 1.5)],
+    "frequency": 5.0,
+    "delay": 0.2,
+    "delay_type": "time",
+    "receiver_locations": spyro.create_transect((-1.45, 0.7), (-1.45, 1.3), 200),
+}
+dictionary["time_axis"] = {
+    "initial_time": 0.0,  # Initial time for event
+    "final_time": final_time,  # Final time for event
+    "dt": 0.001,  # timestep size
+    "amplitude": 1,  # the Ricker has an amplitude of 1.
+    "output_frequency": 100,  # how frequently to output solution to pvds - Perguntar Daiane ''post_processing_frequnecy'
+    "gradient_sampling_frequency": 1,  # how frequently to save solution to RAM    - Perguntar Daiane 'gradient_sampling_frequency'
+}
+dictionary["visualization"] = {
+    "forward_output": False,
+    "forward_output_filename": "results/forward_output.pvd",
+    "fwi_velocity_model_output": False,
+    "velocity_model_filename": None,
+    "gradient_output": False,
+    "gradient_filename": "results/Gradient.pvd",
+    "adjoint_output": False,
+    "adjoint_filename": None,
+    "debug_output": False,
+}
+dictionary["inversion"] = {
+    "perform_fwi": True,  # switch to true to make a FWI
+    "initial_guess_model_file": None,
+    "shot_record_file": None,
+}
+
+def test_real_shot_record_generation_parallel():
+
+    fwi = spyro.FullWaveformInversion(dictionary=dictionary)
+
+    fwi.set_real_mesh(mesh_parameters={"dx": 0.1})
+    center_z = -1.0
+    center_x = 1.0
+    mesh_z = fwi.mesh_z
+    mesh_x = fwi.mesh_x
+    cond = fire.conditional((mesh_z-center_z)**2 + (mesh_x-center_x)**2 < .2**2, 3.0, 2.5)
+
+    fwi.set_real_velocity_model(conditional=cond, output=True, dg_velocity_model=False)
+    fwi.generate_real_shot_record(plot_model=True, save_shot_record=True)
+
+
+def test_realistic_fwi():
+    dictionary["inversion"] = {
+        "real_shot_record_files": "shots/shot_record_",
+    }
+    fwi = spyro.FullWaveformInversion(dictionary=dictionary)
+    fwi.set_guess_mesh(mesh_parameters={"dx": 0.1})
+    fwi.set_guess_velocity_model(constant=2.5)
+    mask_boundaries = {
+        "z_min": -1.3,
+        "z_max": -0.7,
+        "x_min": 0.7,
+        "x_max": 1.3,
+    }
+    fwi.set_gradient_mask(boundaries=mask_boundaries)
+    fwi.run_fwi(vmin=2.5, vmax=3.0, maxiter=5)
+
+
+if __name__ == "__main__":
+    test_real_shot_record_generation_parallel()
+    # test_realistic_fwi()

test_polygon_fwi.py

@@ -5,6 +5,9 @@
 import spyro
 import numpy as np
 import matplotlib.pyplot as plt
+import warnings
+warnings.filterwarnings("ignore")
+
 
 
 def test_real_shot_record_generation_parallel():
@@ -130,9 +133,9 @@ def test_realistic_fwi():
     fwi.set_guess_velocity_model(new_file="velocity_models/case1_sigma10.hdf5")
     mask_boundaries = {
         "z_min": -2.0,
-        "z_max": -0.25,
-        "x_min": 1.0,
-        "x_max": 2.0,
+        "z_max": -0.,
+        "x_min": 0.0,
+        "x_max": 3.0,
     }
     fwi.set_gradient_mask(boundaries=mask_boundaries)
     fwi.run_fwi(vmin=1.5, vmax=3.25, maxiter=5)