improved velocity smoother

spyro/solvers/inversion.py

@@ -353,7 +353,7 @@ def set_real_mesh(
     def set_guess_mesh(
         self,
         user_mesh=None,
-        mesh_parameters=None,
+        mesh_parameters={},
     ):
         """
         Set the mesh for the guess model.

spyro/tools/velocity_smoother.py

@@ -6,7 +6,7 @@
 from SeismicMesh import write_velocity_model
 
 
-def smooth_velocity_field_file(input_filename, output_filename, sigma, show=False, write_hdf5=True):
+def smooth_velocity_field_file(input_filename, output_filename, sigma, show=False, write_hdf5=True, i_limit=None, vp_limit=None, tol=1e-5):
     """Smooths a velocity field using a Gaussian filter.
 
     Parameters
@@ -35,16 +35,24 @@ def smooth_velocity_field_file(input_filename, output_filename, sigma, show=Fals
                 vp[:, index] = trace
     else:
         raise ValueError("Not yet implemented!")
+
+    vp_min = np.min(vp)
+    vp_max = np.max(vp)
+    print(f"Velocity model has minimum vp of {vp_min}, and max of {vp_max}")
 
     vp_smooth = gaussian_filter(vp, sigma)
     ni, nj = np.shape(vp)
+    if i_limit is None:
+        i_limit = 0
+    if vp_limit is None:
+        vp_limit = vp_min
 
     for i in range(ni):
         for j in range(nj):
-            if i < 25:
+            if i < i_limit:
                 vp_smooth[i, j] = vp[i, j]
-            if i < 20:
-                vp_smooth[i, j] = 1.5
+            if vp[i, j] <= vp_limit + tol:
+                vp_smooth[i, j] = vp_min
 
     spec = segyio.spec()
     spec.sorting = 2  # not sure what this means
@@ -73,6 +81,7 @@ def smooth_velocity_field_file(input_filename, output_filename, sigma, show=Fals
         plt.xlabel("x-direction (m)")
         plt.ylabel("z-direction (m)")
         ax.axis("equal")
+        plt.savefig("debug.png")
         plt.show()
 
     if write_hdf5:

test_small_marmousi.py

@@ -74,23 +74,35 @@ def test_real_shot_record_generation_parallel():
 
 
 def test_realistic_fwi():
-    dictionary["inversion"] = {
+    fwi_dictionary = deepcopy(dictionary)
+    fwi_dictionary["inversion"] = {
         "perform_fwi": True,
         "real_shot_record_files": f"shots/shot_record_",
-        "initial_guess_model_file": "velocity_models/initial_guess_15Hz.hdf5",
+        "initial_guess_model_file": "velocity_models/initial_guess_5Hz.hdf5",
     }
+    fwi_dictionary["mesh"]["mesh_type"] = "SeismicMesh"
     fwi = spyro.FullWaveformInversion(dictionary=dictionary)
-    fwi.set_guess_velocity_model(new_file="velocity_models/initial_guess_15Hz.hdf5")
+    # fwi.set_guess_velocity_model(new_file="velocity_models/initial_guess_5Hz.hdf5")
+    fwi.set_guess_mesh()
     mask_boundaries = {
-        "z_min": -1.3,
-        "z_max": -0.7,
-        "x_min": 0.7,
-        "x_max": 1.3,
+        "z_min": -3.5,
+        "z_max": -0.5,
+        "x_min": 3.0,
+        "x_max": 13.0,
     }
     fwi.set_gradient_mask(boundaries=mask_boundaries)
     fwi.run_fwi(vmin=2.5, vmax=3.5, maxiter=60)
 
 
+def test_smoothing_real_model():
+    real_vp_file = "velocity_models/vp_marmousi-ii.segy"
+    sigma = 100
+    output_filename = f"velocity_models/vp_marmousi-ii_sigma{sigma}.segy"
+
+    spyro.tools.velocity_smoother.smooth_velocity_field_file(real_vp_file, output_filename, sigma, show=True, vp_limit=1.5)
+
+
 if __name__ == "__main__":
-    test_real_shot_record_generation_parallel()
+    # test_real_shot_record_generation_parallel()
     # test_realistic_fwi()
+    test_smoothing_real_model()