fwi debugging

cleanup.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 rm *.msh
 rm *.vtk
-rm *.png
+# rm *.png
 rm *.vtu
 rm *.pvtu
 rm *.pvd
@@ -11,8 +11,10 @@ rm *.dat
 rm results/*.vtu
 rm results/*.pvd
 rm results/*.pvtu
-rm shots/*.dat
+# rm shots/*.dat
 rm -rf results/shot*
 rm -rf results/gradient
-rm -rf results/adjoint_shot
-
+rm -rf results/*
+rm -rf control_*/
+rm -rf gradient_*/
+rm -rf initial_velocity_model/

spyro/examples/immersed_polygon.py

@@ -40,8 +40,8 @@
     # every processor) or spatial
 }
 polygon_dictionary["mesh"] = {
-    "Lz": 1.0,  # depth in km - always positive
-    "Lx": 1.0,  # width in km - always positive
+    "Lz": 2.0,  # depth in km - always positive
+    "Lx": 3.0,  # width in km - always positive
     "Ly": 0.0,  # thickness in km - always positive
     "h": 0.05,  # mesh size in km
     "mesh_file": None,
@@ -89,6 +89,7 @@
 }
 polygon_dictionary["polygon_options"] = {
     "water_layer_is_present": True,
+    "water_layer_depth": 0.2,
     "upper_layer": 2.0,
     "middle_layer": 2.5,
     "lower_layer": 3.0,
@@ -131,19 +132,24 @@ def _polygon_velocity_model(self):
         x = self.mesh_x
 
         v0 = 1.5  # water layer
+        water_layer_depth = polygon_dict.get("water_layer_depth", 0.0)
+        water_layer_present = polygon_dict.get("water_layer_is_present", False)
         v1 = polygon_dict["upper_layer"]
         v2 = polygon_dict["middle_layer"]  # background vp (km/s)
         vl = polygon_dict["lower_layer"]  # lower layer (km/s)
         dv = polygon_dict["polygon_layer_perturbation"]*v2  # 30% of perturbation
+        d0 = -water_layer_depth
+        d1 = d0 - 0.14
+        d2 = d1 - 0.2
 
-        if polygon_dict["water_layer_is_present"]:
-            cond = fire.conditional(z >= -0.16, v0, v1)
-            cond = fire.conditional(z <= -0.3, v2, cond)
+        if water_layer_present:
+            cond = fire.conditional(z >= d0, v0, v1)
+            cond = fire.conditional(z <= d1, v2, cond)
         else:
-            cond = fire.conditional(z <= -0.3, v2, v1)
-        cond = fire.conditional(z <= -0.5 - 0.2*x, vl, cond)
+            cond = fire.conditional(z <= d1, v2, v1)
+        cond = fire.conditional(z <= d2 - 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(300*((x-1.5)*(-z-0.7))**2 + ((x-1.5)+(-z-0.7))**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
@@ -196,19 +202,25 @@ def _polygon_velocity_model(self):
         x = self.mesh_x
 
         v0 = 1.5  # water layer
+        water_layer_depth = polygon_dict.get("water_layer_depth", 0.0)
+        water_layer_present = polygon_dict.get("water_layer_is_present", False)
         v1 = polygon_dict["upper_layer"]
         v2 = polygon_dict["middle_layer"]  # background vp (km/s)
         vl = polygon_dict["lower_layer"]  # lower layer (km/s)
         dv = polygon_dict["polygon_layer_perturbation"]*v2  # 30% of perturbation
+        d0 = -water_layer_depth
+        d1 = d0 - 0.14
+        d2 = d1 - 0.2
 
-        if polygon_dict["water_layer_is_present"]:
-            cond = fire.conditional(z >= -0.16, v0, v1)
-            cond = fire.conditional(z <= -0.3, v2, cond)
+        if water_layer_present:
+            cond = fire.conditional(z >= d0, v0, v1)
+            cond = fire.conditional(z <= d1, v2, cond)
         else:
-            cond = fire.conditional(z <= -0.3, v2, v1)
-        cond = fire.conditional(z <= -0.5 - 0.2*x, vl, cond)
+            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
             cond = fire.conditional(z <= middle_of_pad, v0, cond)

spyro/io/basicio.py

@@ -224,10 +224,10 @@ def write_function_to_grid(function, V, grid_spacing):
     x, y = coords.dat.data[:, 0], coords.dat.data[:, 1]
 
     # add buffer to avoid NaN when calling griddata
-    min_x = np.amin(x) + 0.01
-    max_x = np.amax(x) - 0.01
-    min_y = np.amin(y) + 0.01
-    max_y = np.amax(y) - 0.01
+    min_x = np.amin(x) + 0.005
+    max_x = np.amax(x) - 0.005
+    min_y = np.amin(y) + 0.005
+    max_y = np.amax(y) - 0.005
 
     z = function.dat.data[:]
 
@@ -384,7 +384,7 @@ def interpolate(Model, fname, V):
     """
     sd = V.mesh().geometric_dimension()
     m = V.ufl_domain()
-    if Model.abc_active:
+    if Model.abc_active or Model.abc_pad_length > 0.1:
         minz = -Model.length_z - Model.abc_pad_length
         maxz = 0.0
         minx = 0.0 - Model.abc_pad_length

spyro/solvers/backward_time_integration.py

@@ -53,9 +53,9 @@ def backward_wave_propagator_no_pml(Wave_obj, dt=None):
     temp_filename = Wave_obj.forward_output_file
 
     filename, file_extension = temp_filename.split(".")
-    output_filename = "backward." + file_extension
+    # output_filename = "backward." + file_extension
 
-    output = fire.File(output_filename, comm=comm.comm)
+    # output = fire.File(output_filename, comm=comm.comm)
     comm.comm.barrier()
 
     X = fire.Function(Wave_obj.function_space)
@@ -118,8 +118,8 @@ def backward_wave_propagator_no_pml(Wave_obj, dt=None):
                 fire.norm(u_n) < 1
             ), "Numerical instability. Try reducing dt or building the \
                 mesh differently"
-            if Wave_obj.forward_output:
-                output.write(u_n, time=t, name="Pressure")
+            # if Wave_obj.forward_output:
+            #     output.write(u_n, time=t, name="Pressure")
 
             helpers.display_progress(Wave_obj.comm, t)
 

spyro/tools/velocity_smoother.py

@@ -41,8 +41,10 @@ def smooth_velocity_field_file(input_filename, output_filename, sigma, show=Fals
 
     for i in range(ni):
         for j in range(nj):
-            if vp[i, j] < 1.51:
+            if i < 25:
                 vp_smooth[i, j] = vp[i, j]
+            if i < 20:
+                vp_smooth[i, j] = 1.5
 
     spec = segyio.spec()
     spec.sorting = 2  # not sure what this means

test_polygon_fwi.py

@@ -10,22 +10,23 @@
 def test_real_shot_record_generation_parallel():
     dictionary = {}
     dictionary["absorving_boundary_conditions"] = {
-        "pad_length": 2.0,  # True or false
+        "pad_length": 1.0,  # True or false
     }
     dictionary["mesh"] = {
-        "h": 0.05,  # mesh size in km
+        "h": 0.03,  # mesh size in km
     }
     dictionary["polygon_options"] = {
         "water_layer_is_present": True,
+        "water_layer_depth": 0.2,
         "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), 4),
+        "source_locations": spyro.create_transect((-0.1, 1.0), (-0.1, 2.0), 7),
         "frequency": 5.0,
-        "receiver_locations": spyro.create_transect((-0.16, 0.1), (-0.16, 0.9), 100),
+        "receiver_locations": spyro.create_transect((-0.15, 1.0), (-0.15, 2.0), 100),
     }
     dictionary["visualization"] = {
         "debug_output": True,
@@ -41,39 +42,42 @@ def test_real_shot_record_generation_parallel():
     fwi = spyro.examples.Polygon_acoustic_FWI(dictionary=dictionary, periodic=True)
     fwi.generate_real_shot_record(plot_model=True, save_shot_record=True)
 
-    dictionary["inversion"] = {
-        "real_shot_record_files": "shots/shot_record_",
-    }
-    fwi2 = spyro.examples.Polygon_acoustic_FWI(dictionary=dictionary, periodic=True)
+    # dictionary["inversion"] = {
+    #     "real_shot_record_files": "shots/shot_record_",
+    # }
+    # fwi2 = spyro.examples.Polygon_acoustic_FWI(dictionary=dictionary, periodic=True)
 
-    max_value = np.max(fwi2.real_shot_record[:, fwi2.comm.ensemble_comm.rank*33])
-    test_core = np.isclose(max_value, 0.184, atol=1e-2)
+    # max_value = np.max(fwi2.real_shot_record[:, fwi2.comm.ensemble_comm.rank*33])
+    # test_core = np.isclose(max_value, 0.184, atol=1e-2)
 
-    test_core_all = fwi2.comm.ensemble_comm.allgather(test_core)
-    test = all(test_core_all)
+    # test_core_all = fwi2.comm.ensemble_comm.allgather(test_core)
+    # test = all(test_core_all)
 
-    print(f"Correctly loaded shots: {test}", flush=True)
+    # print(f"Correctly loaded shots: {test}", flush=True)
 
-    assert test
+    # assert test
 
 
 def test_velocity_smoother_in_fwi():
     dictionary = {}
     dictionary["absorving_boundary_conditions"] = {
-        "pad_length": 2.0,  # True or false
+        "pad_length": 1.0,  # True or false
     }
     dictionary["mesh"] = {
-        "h": 0.01,  # mesh size in km
+        "h": 0.03,  # mesh size in km
     }
     dictionary["polygon_options"] = {
         "water_layer_is_present": True,
+        "water_layer_depth": 0.2,
         "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),
+        "source_locations": spyro.create_transect((-0.1, 1.0), (-0.1, 2.0), 1),
+        "frequency": 5.0,
+        "receiver_locations": spyro.create_transect((-0.15, 1.0), (-0.15, 2.0), 100),
     }
     fwi = spyro.examples.Polygon_acoustic_FWI(dictionary=dictionary, periodic=True)
     spyro.io.create_segy(
@@ -87,6 +91,54 @@ def test_velocity_smoother_in_fwi():
     plt.close()
 
 
+def test_realistic_fwi():
+    dictionary = {}
+    dictionary["absorving_boundary_conditions"] = {
+        "pad_length": 1.0,  # True or false
+    }
+    dictionary["mesh"] = {
+        "h": 0.05,  # mesh size in km
+    }
+    dictionary["polygon_options"] = {
+        "water_layer_is_present": True,
+        "water_layer_depth": 0.2,
+        "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, 1.0), (-0.1, 2.0), 7),
+        "frequency": 5.0,
+        "receiver_locations": spyro.create_transect((-0.15, 1.0), (-0.15, 2.0), 100),
+    }
+    dictionary["visualization"] = {
+        "debug_output": True,
+    }
+    dictionary["time_axis"] = {
+        "final_time": 1.0,  # Final time for event
+        "dt": 0.0005,  # 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"] = {
+        "real_shot_record_files": "shots/shot_record_",
+    }
+    fwi = spyro.examples.Polygon_acoustic_FWI(dictionary=dictionary, periodic=True)
+    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,
+    }
+    fwi.set_gradient_mask(boundaries=mask_boundaries)
+    fwi.run_fwi(vmin=1.5, vmax=3.25, maxiter=5)
+
+
 if __name__ == "__main__":
     # test_real_shot_record_generation_parallel()
-    test_velocity_smoother_in_fwi()
+    # test_velocity_smoother_in_fwi()
+    test_realistic_fwi()