started testing

case1.py

@@ -0,0 +1,61 @@
+# from mpi4py.MPI import COMM_WORLD
+# import debugpy
+# debugpy.listen(3000 + COMM_WORLD.rank)
+# debugpy.wait_for_client()
+import spyro
+import matplotlib.pyplot as plt
+dictionary = {}
+dictionary["absorving_boundary_conditions"] = {
+    "pad_length": 2.0,  # True or false
+}
+dictionary["mesh"] = {
+    "h": 0.05,  # mesh size in km
+}
+dictionary["polygon_options"] = {
+    "water_layer_is_present": True,
+    "upper_layer": 2.0,
+    "middle_layer": 2.5,
+    "lower_layer": 3.0,
+    "polygon_layer_perturbation": 0.3,
+}
+dictionary["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": spyro.create_transect((-0.1, 0.1), (-0.1, 0.9), 1),
+    "frequency": 5.0,
+    "delay": 1.0,
+    "receiver_locations": spyro.create_transect((-0.16, 0.1), (-0.16, 0.9), 100),
+}
+dictionary["visualization"] = {
+    "debug_output": True,
+}
+dictionary["time_axis"] = {
+    "initial_time": 0.0,  # Initial time for event
+    "final_time": 1.0,  # Final time for event
+    "dt": 0.0001,  # timestep size
+    "amplitude": 1,  # the Ricker has an amplitude of 1.
+    "output_frequency": 500,  # how frequently to output solution to pvds
+    # how frequently to save solution to RAM
+    "gradient_sampling_frequency": 1,
+}
+dictionary["inversion"] = {
+    "initial_guess_model_file": "velocity_models/case1_sigma10.hdf5",
+}
+fwi = spyro.examples.Polygon_acoustic_FWI(dictionary=dictionary, periodic=True)
+fwi.generate_real_shot_record(plot_model=True)
+fwi.initial_velocity_model = None
+fwi.set_guess_velocity_model(new_file="velocity_models/case1_sigma10.hdf5")
+# spyro.io.save_shots(fwi)
+# spyro.io.load_shots(fwi)
+# fwi.real_shot_record = fwi.forward_solution_receivers
+fwi.run_fwi(vmin=1.5, vmax=3.25, maxiter=5)
+# spyro.io.create_segy(
+#     fwi.initial_velocity_model,
+#     fwi.function_space,
+#     10.0/1000.0,
+#     "velocity_models/true_case1.segy",
+#     )
+
+# plt.close()
+# spyro.tools.velocity_smoother.smooth_velocity_field_file("velocity_models/true_case1.segy", "velocity_models/case1_sigma10.segy", 10, show=True)
+# plt.savefig("velocity_models/case1_sigma10.png")
+print("END")

case1_v2.py

@@ -0,0 +1,47 @@
+# from mpi4py.MPI import COMM_WORLD
+# import debugpy
+# debugpy.listen(3000 + COMM_WORLD.rank)
+# debugpy.wait_for_client()
+import spyro
+import matplotlib.pyplot as plt
+dictionary = {}
+dictionary["absorving_boundary_conditions"] = {
+    "pad_length": 2.0,  # True or false
+}
+dictionary["mesh"] = {
+    "h": 0.05,  # mesh size in km
+}
+dictionary["polygon_options"] = {
+    "water_layer_is_present": True,
+    "upper_layer": 2.0,
+    "middle_layer": 2.5,
+    "lower_layer": 3.0,
+    "polygon_layer_perturbation": 0.3,
+}
+dictionary["acquisition"] = {
+    "source_type": "ricker",
+    "source_locations": spyro.create_transect((-0.1, 0.1), (-0.1, 0.9), 1),
+    "frequency": 5.0,
+    "delay": 1.0,
+    "receiver_locations": spyro.create_transect((-0.16, 0.1), (-0.16, 0.9), 100),
+}
+dictionary["visualization"] = {
+    "debug_output": True,
+}
+dictionary["time_axis"] = {
+    "initial_time": 0.0,  # Initial time for event
+    "final_time": 1.0,  # Final time for event
+    "dt": 0.0001,  # timestep size
+    "amplitude": 1,  # the Ricker has an amplitude of 1.
+    "output_frequency": 500,  # how frequently to output solution to pvds
+    # how frequently to save solution to RAM
+    "gradient_sampling_frequency": 1,
+}
+dictionary["inversion"] = {
+    "initial_guess_model_file": "velocity_models/case1_sigma10.hdf5",
+}
+fwi = spyro.examples.Polygon_acoustic_FWI(dictionary=dictionary, periodic=True)
+fwi.generate_real_shot_record(plot_model=True)
+spyro.io.save_shots(fwi)
+
+print("END")

spyro/examples/immersed_polygon.py

@@ -144,7 +144,10 @@ def _polygon_velocity_model(self):
         cond = fire.conditional(z <= -0.5 - 0.2*x, vl, 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)
-        cond = fire.conditional(z <= -2.0,v0 , cond)
+
+        if self.abc_pad_length > 0.0:
+            middle_of_pad = -self.length_z - self.abc_pad_length*0.5
+            cond = fire.conditional(z <= middle_of_pad, v0, cond)
 
         self.set_initial_velocity_model(conditional=cond, dg_velocity_model=False)
         return None
@@ -206,7 +209,9 @@ def _polygon_velocity_model(self):
         cond = fire.conditional(z <= -0.5 - 0.2*x, vl, 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)
-        cond = fire.conditional(z <= -2.0,v0 , cond)
+        if self.abc_pad_length > 0.0:
+            middle_of_pad = -self.length_z - self.abc_pad_length*0.5
+            cond = fire.conditional(z <= middle_of_pad, v0, cond)
 
         self.set_initial_velocity_model(conditional=cond, dg_velocity_model=False)
         return None

spyro/io/basicio.py

@@ -307,7 +307,7 @@ def rebuild_empty_forward_solution(wave, time_steps):
 
 
 @ensemble_save_or_load
-def load_shots(Wave_obj, file_name=None, shot_ids=0):
+def load_shots(Wave_obj, file_name="shots/shot_record_", shot_ids=0):
     """Load a `pickle` to a `numpy.ndarray`.
 
     Parameters

spyro/io/boundary_layer_io.py

@@ -49,7 +49,7 @@ def __init__(self, abc_dictionary):
             self.abc_exponent = None
             self.abc_cmax = None
             self.abc_R = None
-            self.abc_pad_length = 0.0
+            self.abc_pad_length = self.dictionary.get("pad_length", 0.0)
             self.abc_boundary_layer_type = None
             pass
         elif self.dictionary["damping_type"] == "PML":

test_polygon_example.py

@@ -0,0 +1,96 @@
+import spyro
+import numpy as np
+
+
+def test_polygon_vp():
+    dictionary = {}
+    dictionary["polygon_options"] = {
+        "water_layer_is_present": True,
+        "upper_layer": 2.0,
+        "middle_layer": 2.5,
+        "lower_layer": 3.0,
+        "polygon_layer_perturbation": 0.3,
+    }
+    dictionary["absorving_boundary_conditions"] = {
+        "pad_length": 1.0,  # True or false
+    }
+    wave = spyro.examples.Polygon_acoustic(dictionary=dictionary, periodic=True)
+
+    # Check water layer velocity
+    test_locations = [
+        (-0.1, 0.5),  # Water layer p1
+        (-0.05, 0.7),  # Water layer p2
+        (-0.22, 0.1),  # Upper layer p1
+        (-0.25, 0.6), # Upper layer p2
+        (-0.50, 0.1), # Middle layer p1
+        (-0.55, 0.1), # Bottom layer p1
+        (-0.57, 0.2), # Bottom layer p2
+        (-0.3, 0.5),  # polygon p1
+        (-0.4, 0.6),  # polygon p2
+        (-1.3, 0.5),  # pad before change
+        (-1.6, 0.5),  # pad after change
+    ]
+    expected_values = [
+        1.5,
+        1.5,
+        dictionary["polygon_options"]["upper_layer"],
+        dictionary["polygon_options"]["upper_layer"],
+        dictionary["polygon_options"]["middle_layer"],
+        dictionary["polygon_options"]["lower_layer"],
+        dictionary["polygon_options"]["lower_layer"],
+        dictionary["polygon_options"]["middle_layer"]*(1+dictionary["polygon_options"]["polygon_layer_perturbation"]),
+        dictionary["polygon_options"]["middle_layer"]*(1+dictionary["polygon_options"]["polygon_layer_perturbation"]),
+        dictionary["polygon_options"]["lower_layer"],
+        1.5,
+    ]
+
+    # Check upper layer
+    test_array = np.isclose(wave.initial_velocity_model.at(test_locations), expected_values)
+    test = test_array.all()
+
+    print(f"All points arrive at expected values: {test}", flush=True)
+    assert test
+
+
+def test_real_shot_record_generation_for_polygon():
+    dictionary = {}
+    dictionary["absorving_boundary_conditions"] = {
+        "pad_length": 2.0,  # True or false
+    }
+    dictionary["mesh"] = {
+        "h": 0.05,  # mesh size in km
+    }
+    dictionary["polygon_options"] = {
+        "water_layer_is_present": True,
+        "upper_layer": 2.0,
+        "middle_layer": 2.5,
+        "lower_layer": 3.0,
+        "polygon_layer_perturbation": 0.3,
+    }
+    dictionary["acquisition"] = {
+        "source_locations": spyro.create_transect((-0.1, 0.1), (-0.1, 0.9), 1),
+        "frequency": 5.0,
+        "receiver_locations": spyro.create_transect((-0.16, 0.1), (-0.16, 0.9), 100),
+    }
+    dictionary["visualization"] = {
+        "debug_output": True,
+    }
+    dictionary["time_axis"] = {
+        "final_time": 1.0,  # Final time for event
+        "dt": 0.0001,  # timestep size
+        "amplitude": 1,  # the Ricker has an amplitude of 1.
+        "output_frequency": 500,  # how frequently to output solution to pvds
+        # how frequently to save solution to RAM
+        "gradient_sampling_frequency": 1,
+    }
+    fwi = spyro.examples.Polygon_acoustic_FWI(dictionary=dictionary, periodic=True)
+    fwi.generate_real_shot_record(plot_model=True)
+    spyro.io.save_shots(fwi)
+
+    print("END")
+
+
+
+if __name__ == "__main__":
+    # test_polygon_vp()
+    test_real_shot_record_generation_for_polygon()