From 63b7c278cbf551d1140b062bfde82c080d29f40d Mon Sep 17 00:00:00 2001 From: olender Date: Fri, 22 Nov 2024 18:08:38 -0300 Subject: [PATCH] linting and debugging --- save_new_shot_records.py | 2 +- shot_filters.py | 27 ++-- spyro/examples/camembert.py | 1 - spyro/examples/immersed_polygon.py | 3 +- spyro/examples/rectangle.py | 1 + spyro/io/basicio.py | 3 +- spyro/io/model_parameters.py | 10 +- spyro/meshing/meshing_functions.py | 4 +- spyro/solvers/inversion.py | 8 +- .../time_integration_central_difference.py | 2 +- spyro/solvers/wave.py | 12 +- spyro/tools/velocity_smoother.py | 2 +- spyro/utils/utils.py | 2 - temp_test_fwi.py | 140 ++++++++++++++++++ test/test_serialshot_fwi.py | 2 +- test_camembert_fwi.py | 7 +- test_parallel/test_gradient_serialshots.py | 2 +- test_polygon_example.py | 8 +- test_polygon_fwi.py | 1 - test_queijo_minas.py | 15 +- test_small_marmousi.py | 2 +- 21 files changed, 199 insertions(+), 55 deletions(-) create mode 100644 temp_test_fwi.py diff --git a/save_new_shot_records.py b/save_new_shot_records.py index 1e9ddc43..f6dfff2c 100644 --- a/save_new_shot_records.py +++ b/save_new_shot_records.py @@ -21,4 +21,4 @@ print("PAUSE") with open(file_name, "wb") as f: - pickle.dump(sliced_array, f) \ No newline at end of file + pickle.dump(sliced_array, f) diff --git a/shot_filters.py b/shot_filters.py index 81776023..5d14282a 100644 --- a/shot_filters.py +++ b/shot_filters.py @@ -1,6 +1,6 @@ import spyro import numpy as np -from scipy.signal import butter, filtfilt, sosfilt +from scipy.signal import butter, filtfilt, sosfilt from scipy.signal import sosfilt import sys @@ -8,30 +8,31 @@ filter_frequency = 7.0 -def filter_shot(shot, cutoff, fs, filter_type = 'butter'): + +def filter_shot(shot, cutoff, fs, filter_type='butter'): if filter_type == 'butter': - return butter_filter(shot,cutoff, fs) + return butter_filter(shot, cutoff, fs) def butter_filter(shot, cutoff, fs, order=1): - """ Low-pass filter the shot record with sampling-rate fs Hz and cutoff freq. Hz """ - - nyq = 0.5*fs # Nyquist Frequency + + nyq = 0.5*fs # Nyquist Frequency normal_cutoff = (cutoff) / nyq - - # Get the filter coefficients + + # Get the filter coefficients b, a = butter(order, normal_cutoff, btype="low", analog=False) - + nc, nr = np.shape(shot) for rec in range(nr): - shot[:,rec] = filtfilt(b, a, shot[:,rec]) - + shot[:, rec] = filtfilt(b, a, shot[:, rec]) + return shot + frequency = 7.0 dt = 0.0001 degree = 4 @@ -91,7 +92,7 @@ def butter_filter(shot, cutoff, fs, order=1): spyro.io.load_shots(fwi, file_name="shots/shot_record_") shots = fwi.forward_solution_receivers shots *= 5.455538535049624 -print(f'Applying {filter_type} filter for {filter_frequency}Hz', flush = True) -p_filtered = filter_shot(shots, filter_frequency, fs, filter_type = filter_type) +print(f'Applying {filter_type} filter for {filter_frequency}Hz', flush=True) +p_filtered = filter_shot(shots, filter_frequency, fs, filter_type=filter_type) shot_filename = f"shots/shot_record_{filter_frequency}_" spyro.io.save_shots(fwi, file_name=shot_filename) diff --git a/spyro/examples/camembert.py b/spyro/examples/camembert.py index 7169ee76..7557e0dd 100644 --- a/spyro/examples/camembert.py +++ b/spyro/examples/camembert.py @@ -211,4 +211,3 @@ def _camembert_velocity_model(self): ) self.set_initial_velocity_model(conditional=cond, dg_velocity_model=False) return None - diff --git a/spyro/examples/immersed_polygon.py b/spyro/examples/immersed_polygon.py index b4ea7e28..03a419ec 100644 --- a/spyro/examples/immersed_polygon.py +++ b/spyro/examples/immersed_polygon.py @@ -165,6 +165,7 @@ def _polygon_velocity_model(self): "shot_record_file": None, } + class Polygon_acoustic_FWI(Rectangle_acoustic_FWI): """polygon model. This class is a child of the Example_model class. @@ -226,4 +227,4 @@ def _polygon_velocity_model(self): cond = fire.conditional(z <= middle_of_pad, v0, cond) self.set_initial_velocity_model(conditional=cond, dg_velocity_model=False) - return None \ No newline at end of file + return None diff --git a/spyro/examples/rectangle.py b/spyro/examples/rectangle.py index 61e26199..89cd3600 100644 --- a/spyro/examples/rectangle.py +++ b/spyro/examples/rectangle.py @@ -208,6 +208,7 @@ class Rectangle_acoustic_FWI(Example_model_acoustic_FWI): If True, the mesh will be periodic in all directions. The default is False. """ + def __init__( self, dictionary=None, diff --git a/spyro/io/basicio.py b/spyro/io/basicio.py index dd7a7467..32157cd3 100644 --- a/spyro/io/basicio.py +++ b/spyro/io/basicio.py @@ -154,7 +154,6 @@ def wrapper(*args, **kwargs): J_total[0] /= comm.comm.size elif args[0].parallelism_type == "spatial" and args[0].number_of_sources > 1: - num = args[0].number_of_sources residual_list = args[1] J_total = np.zeros((1)) @@ -205,7 +204,7 @@ def wrapper(*args, **kwargs): kwargs, misfit=current_misfit, ) - ) + ) grad_total += grad grad_total /= num diff --git a/spyro/io/model_parameters.py b/spyro/io/model_parameters.py index a2d2ec48..d88daee4 100644 --- a/spyro/io/model_parameters.py +++ b/spyro/io/model_parameters.py @@ -1,11 +1,12 @@ import numpy as np import uuid -from mpi4py import MPI -from firedrake import COMM_WORLD +from mpi4py import MPI # noqa:F401 +from firedrake import COMM_WORLD # noqa: import warnings from .. import io from .. import utils from .. import meshing +from ..meshing.meshing_functions import cells_per_wavelength # default_optimization_parameters = { # "General": {"Secant": {"Type": "Limited-Memory BFGS", @@ -518,7 +519,6 @@ def _sanitize_comm(self, comm): if self.parallelism_type == "custom": self.shot_ids_per_propagation = dictionary["parallelism"]["shot_ids_per_propagation"] elif self.parallelism_type == "automatic": - available_cores = COMM_WORLD.size self.shot_ids_per_propagation = [[i] for i in range(0, self.number_of_sources)] elif self.parallelism_type == "spatial": self.shot_ids_per_propagation = [[i] for i in range(0, self.number_of_sources)] @@ -606,7 +606,7 @@ def _sanitize_optimization_and_velocity_for_fwi(self): self.initial_velocity_model_file = dictionary["inversion"][ "initial_guess_model_file" ] - except: + except KeyError: self.initial_velocity_model_file = None self.fwi_output_folder = "fwi/" self.control_output_file = self.fwi_output_folder + "control" @@ -638,6 +638,7 @@ def _sanitize_optimization_and_velocity_for_fwi(self): if "shot_record_file" in dictionary["inversion"]: if dictionary["inversion"]["shot_record_file"] is not None: self.real_shot_record = np.load(dictionary["inversion"]["shot_record_file"]) + def _sanitize_optimization_and_velocity_without_fwi(self): dictionary = self.input_dictionary if "synthetic_data" in dictionary: @@ -721,6 +722,7 @@ def set_mesh( mesh_parameters.setdefault("degree", self.degree) mesh_parameters.setdefault("velocity_model_file", self.initial_velocity_model_file) mesh_parameters.setdefault("cell_type", self.cell_type) + print(f"Method: {self.method}, Degree: {self.degree}, Dimension: {self.dimension}") mesh_parameters.setdefault("cells_per_wavelength", cells_per_wavelength(self.method, self.degree, self.dimension)) self._set_mesh_length( diff --git a/spyro/meshing/meshing_functions.py b/spyro/meshing/meshing_functions.py index fb2bba8a..05b8b0b8 100644 --- a/spyro/meshing/meshing_functions.py +++ b/spyro/meshing/meshing_functions.py @@ -17,9 +17,9 @@ def cells_per_wavelength(method, degree, dimension): 'mlt3tet': 3.72, } - if dimension == 2 and (method == 'MLT' or method == 'CG'): + if dimension == 2 and (method == 'mass_lumped_triangle'): cell_type = 'tri' - if dimension == 3 and (method == 'MLT' or method == 'CG'): + if dimension == 3 and (method == 'mass_lumped_triangle'): cell_type = 'tet' key = method.lower()+str(degree)+cell_type diff --git a/spyro/solvers/inversion.py b/spyro/solvers/inversion.py index 8eb4a4df..4f7d6d60 100644 --- a/spyro/solvers/inversion.py +++ b/spyro/solvers/inversion.py @@ -3,14 +3,12 @@ from scipy.optimize import minimize as scipy_minimize from mpi4py import MPI # noqa: F401 import numpy as np -from copy import deepcopy import resource from .acoustic_wave import AcousticWave from ..utils import compute_functional from ..utils import Gradient_mask_for_pml, Mask from ..plots import plot_model as spyro_plot_model -from ..io.basicio import ensemble_shot_record from ..io.basicio import switch_serial_shot from ..io.basicio import load_shots, save_shots @@ -200,7 +198,7 @@ def calculate_misfit(self, c=None): if self.parallelism_type == "spatial" and self.number_of_sources > 1: misfit_list = [] guess_shot_record_list = [] - for snum in range (self.number_of_sources): + for snum in range(self.number_of_sources): switch_serial_shot(self, snum) guess_shot_record_list.append(self.forward_solution_receivers) misfit_list.append(self.real_shot_record[snum] - self.forward_solution_receivers) @@ -396,7 +394,7 @@ def get_functional(self, c=None): print(f"Functional: {Jm} at iteration: {self.current_iteration}", flush=True) with open("functional_values.txt", "a") as file: file.write(f"Iteration: {self.current_iteration}, Functional: {Jm}\n") - + with open("peak_memory.txt", "a") as file: file.write(f"Peak memory usage: {peak_memory_mb:.2f} MB \n") @@ -619,7 +617,7 @@ def forward_solve(self): super().forward_solve() if self.parallelism_type == "spatial" and self.number_of_sources > 1: real_shot_record_list = [] - for snum in range (self.number_of_sources): + for snum in range(self.number_of_sources): switch_serial_shot(self, snum) real_shot_record_list.append(self.receivers_output) self.real_shot_record = real_shot_record_list diff --git a/spyro/solvers/time_integration_central_difference.py b/spyro/solvers/time_integration_central_difference.py index d6880cc5..882956d4 100644 --- a/spyro/solvers/time_integration_central_difference.py +++ b/spyro/solvers/time_integration_central_difference.py @@ -24,7 +24,7 @@ def central_difference(wave, source_ids=[0]): wave.sources.current_sources = source_ids rhs_forcing = fire.Cofunction(wave.function_space.dual()) - wave.field_logger.start_logging(source_id) + wave.field_logger.start_logging(source_ids) wave.comm.comm.barrier() diff --git a/spyro/solvers/wave.py b/spyro/solvers/wave.py index 39e449ba..3f8d4512 100644 --- a/spyro/solvers/wave.py +++ b/spyro/solvers/wave.py @@ -118,9 +118,9 @@ def set_mesh( self, user_mesh=None, mesh_parameters={}, - ): - """ - Set the mesh for the solver. + ): + """ + Set the mesh for the solver. Args: user_mesh (optional): User-defined mesh. Defaults to None. @@ -373,7 +373,7 @@ def update_source_expression(self, t): pass @ensemble_propagator - def wave_propagator(self, dt=None, final_time=None, source_num=0): + def wave_propagator(self, dt=None, final_time=None, source_nums=[0]): """Propagates the wave forward in time. Currently uses central differences. @@ -398,8 +398,8 @@ def wave_propagator(self, dt=None, final_time=None, source_num=0): if dt is not None: self.dt = dt - self.current_source = source_num - usol, usol_recv = time_integrator(self, source_num) + self.current_sources = source_nums + usol, usol_recv = time_integrator(self, source_nums) return usol, usol_recv diff --git a/spyro/tools/velocity_smoother.py b/spyro/tools/velocity_smoother.py index 164d824b..5489082a 100644 --- a/spyro/tools/velocity_smoother.py +++ b/spyro/tools/velocity_smoother.py @@ -35,7 +35,7 @@ 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}") diff --git a/spyro/utils/utils.py b/spyro/utils/utils.py index 998373f1..aebc5d22 100644 --- a/spyro/utils/utils.py +++ b/spyro/utils/utils.py @@ -46,7 +46,6 @@ def compute_functional(Wave_object, residual): """ num_receivers = Wave_object.number_of_receivers dt = Wave_object.dt - comm = Wave_object.comm J = 0 for rn in range(num_receivers): @@ -95,7 +94,6 @@ def mpi_init(model): num_cores_per_propagation = available_cores elif model.parallelism_type == "custom": shot_ids_per_propagation = model.shot_ids_per_propagation - num_max_shots_per_core = max(len(sublist) for sublist in shot_ids_per_propagation) num_propagations = len(shot_ids_per_propagation) num_cores_per_propagation = available_cores / num_propagations diff --git a/temp_test_fwi.py b/temp_test_fwi.py new file mode 100644 index 00000000..f844f56b --- /dev/null +++ b/temp_test_fwi.py @@ -0,0 +1,140 @@ +import numpy as np +import firedrake as fire +import spyro +import pytest + + +def is_rol_installed(): + try: + import ROL + return True + except ImportError: + return False + + +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), 1), + # "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_fwi(load_real_shot=False, use_rol=False): + """ + Run the Full Waveform Inversion (FWI) test. + + Parameters + ---------- + load_real_shot (bool, optional): Whether to load a real shot record or not. Defaults to False. + """ + + # Setting up to run synthetic real problem + if load_real_shot is False: + FWI_obj = spyro.FullWaveformInversion(dictionary=dictionary) + + FWI_obj.set_real_mesh(mesh_parameters={"dx": 0.1}) + center_z = -1.0 + center_x = 1.0 + mesh_z = FWI_obj.mesh_z + mesh_x = FWI_obj.mesh_x + cond = fire.conditional((mesh_z-center_z)**2 + (mesh_x-center_x)**2 < .2**2, 3.0, 2.5) + + FWI_obj.set_real_velocity_model(conditional=cond, output=True, dg_velocity_model=False) + FWI_obj.generate_real_shot_record( + plot_model=True, + model_filename="True_experiment.png", + abc_points=[(-0.5, 0.5), (-1.5, 0.5), (-1.5, 1.5), (-0.5, 1.5)] + ) + np.save("real_shot_record", FWI_obj.real_shot_record) + + else: + dictionary["inversion"]["shot_record_file"] = "real_shot_record.npy" + FWI_obj = spyro.FullWaveformInversion(dictionary=dictionary) + + # Setting up initial guess problem + FWI_obj.set_guess_mesh(mesh_parameters={"dx": 0.1}) + FWI_obj.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_obj.set_gradient_mask(boundaries=mask_boundaries) + if use_rol: + FWI_obj.run_fwi_rol(vmin=2.5, vmax=3.0, maxiter=2) + else: + FWI_obj.run_fwi(vmin=2.5, vmax=3.0, maxiter=5) + + # simple mask test + grad_test = FWI_obj.gradient + test0 = np.isclose(grad_test.at((-0.1, 0.1)), 0.0) + print(f"PML looks masked: {test0}", flush=True) + test1 = np.abs(grad_test.at((-1.0, 1.0))) > 1e-5 + print(f"Center looks unmasked: {test1}", flush=True) + + # quick look at functional and if it reduced + test2 = FWI_obj.functional < 1e-3 + print(f"Last functional small: {test2}", flush=True) + test3 = FWI_obj.functional_history[-1]/FWI_obj.functional_history[0] < 1e-2 + print(f"Considerable functional reduction during test: {test3}", flush=True) + + print("END", flush=True) + assert all([test0, test1, test2, test3]) + + +@pytest.mark.skipif(not is_rol_installed(), reason="ROL is not installed") +def test_fwi_with_rol(load_real_shot=False, use_rol=True): + test_fwi(load_real_shot=load_real_shot, use_rol=use_rol) + + +if __name__ == "__main__": + test_fwi(load_real_shot=False) + test_fwi_with_rol() diff --git a/test/test_serialshot_fwi.py b/test/test_serialshot_fwi.py index dfece0f5..10c883a2 100644 --- a/test/test_serialshot_fwi.py +++ b/test/test_serialshot_fwi.py @@ -10,7 +10,7 @@ def is_rol_installed(): try: - import ROL + import ROL # noqa:F401 return True except ImportError: return False diff --git a/test_camembert_fwi.py b/test_camembert_fwi.py index b4cc0cbe..3dadda11 100644 --- a/test_camembert_fwi.py +++ b/test_camembert_fwi.py @@ -10,6 +10,7 @@ frequency = float(sys.argv[1]) degree = int(sys.argv[2]) + def cells_per_wavelength(degree): cell_per_wavelength_dictionary = { 'kmv2tri': 7.20, @@ -24,9 +25,10 @@ def cells_per_wavelength(degree): cell_type = 'tri' key = 'kmv'+str(degree)+cell_type - + return cell_per_wavelength_dictionary.get(key) + cpw = cells_per_wavelength(degree) final_time = 0.9 # dx = 2.5/(frequency*cpw) @@ -83,8 +85,9 @@ def cells_per_wavelength(degree): "shot_record_file": None, } + def test_real_shot_record_generation_parallel(): - + fwi = spyro.FullWaveformInversion(dictionary=dictionary) fwi.set_real_mesh(mesh_parameters={"dx": dx}) diff --git a/test_parallel/test_gradient_serialshots.py b/test_parallel/test_gradient_serialshots.py index 3a8775db..54436bf9 100644 --- a/test_parallel/test_gradient_serialshots.py +++ b/test_parallel/test_gradient_serialshots.py @@ -91,7 +91,7 @@ def get_gradient(parallelism_type, points): spyro.io.switch_serial_shot(Wave_obj_exact, source_id) spyro.io.switch_serial_shot(Wave_obj_guess, source_id) misfit_list.append(Wave_obj_exact.forward_solution_receivers - Wave_obj_guess.forward_solution_receivers) - misfit= misfit_list + misfit = misfit_list gradient = Wave_obj_guess.gradient_solve(misfit=misfit) Wave_obj_guess.comm.comm.barrier() diff --git a/test_polygon_example.py b/test_polygon_example.py index 401e17bb..c09b78f8 100644 --- a/test_polygon_example.py +++ b/test_polygon_example.py @@ -21,10 +21,10 @@ def test_polygon_vp(): (-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.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 diff --git a/test_polygon_fwi.py b/test_polygon_fwi.py index bd421e1c..f6b71d73 100644 --- a/test_polygon_fwi.py +++ b/test_polygon_fwi.py @@ -9,7 +9,6 @@ warnings.filterwarnings("ignore") - def test_real_shot_record_generation_parallel(): dictionary = {} dictionary["absorving_boundary_conditions"] = { diff --git a/test_queijo_minas.py b/test_queijo_minas.py index c5efe3c9..ac45422c 100644 --- a/test_queijo_minas.py +++ b/test_queijo_minas.py @@ -10,6 +10,7 @@ degree = int(sys.argv[2]) frequency = float(sys.argv[1]) + def cells_per_wavelength(degree): cell_per_wavelength_dictionary = { 'kmv2tri': 7.20, @@ -24,9 +25,10 @@ def cells_per_wavelength(degree): cell_type = 'tri' key = 'kmv'+str(degree)+cell_type - + return cell_per_wavelength_dictionary.get(key) + cpw = cells_per_wavelength(degree) final_time = 0.9 # dx = 2.5/(frequency*cpw) @@ -84,8 +86,9 @@ def cells_per_wavelength(degree): "shot_record_file": None, } + def test_real_shot_record_generation_parallel(): - + fwi = spyro.FullWaveformInversion(dictionary=dictionary) fwi.set_real_mesh(mesh_parameters={"dx": dx}) @@ -93,12 +96,12 @@ def test_real_shot_record_generation_parallel(): center_x = 1.0 mesh_z = fwi.mesh_z mesh_x = fwi.mesh_x - square_top_z = -0.9 - square_bot_z = -1.1 - square_left_x = 0.9 + square_top_z = -0.9 + square_bot_z = -1.1 + square_left_x = 0.9 square_right_x = 1.1 cond = fire.conditional((mesh_z-center_z)**2 + (mesh_x-center_x)**2 < .2**2, 3.0, 2.5) - cond = fire.conditional( + cond = fire.conditional( fire.And( fire.And(mesh_z < square_top_z, mesh_z > square_bot_z), fire.And(mesh_x > square_left_x, mesh_x < square_right_x) diff --git a/test_small_marmousi.py b/test_small_marmousi.py index 49dd1efc..f34d4843 100644 --- a/test_small_marmousi.py +++ b/test_small_marmousi.py @@ -66,7 +66,7 @@ def test_real_shot_record_generation_parallel(): real_dictionary = deepcopy(dictionary) real_dictionary["mesh"]["mesh_file"] = "meshes/real5hz.msh" - + real_wave = spyro.AcousticWave(dictionary=real_dictionary) real_wave.set_initial_velocity_model(new_file="velocity_models/vp_marmousi-ii.hdf5") real_wave.forward_solve()