Skip to content

Commit

Permalink
debugging
Browse files Browse the repository at this point in the history
  • Loading branch information
Olender committed Aug 23, 2024
1 parent 2b2d6f8 commit 9c93b49
Show file tree
Hide file tree
Showing 6 changed files with 167 additions and 10 deletions.
39 changes: 39 additions & 0 deletions full_fwi.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=40
#SBATCH --partition=amd_large_2
#SBATCH --time=5:00:00
#SBATCH --job-name=MS_marmousi
#SBATCH --output=%x.%j.out
#SBATCH --error=%x.%j.err
#SBATCH --exclusive

echo -e "\n## Job started at $(date +'%d-%m-%Y as %T') #####################\n"
echo -e "\n## Jobs activated by $USER: \n"
squeue -a --user=$USER
echo -e "\n## Execution node: $(hostname -s) \n"
echo -e "\n## Number of tasks per job: $SLURM_NTASKS \n"
#########################################
##------- Start of job ----- #
#########################################
## Configure the execution environment


##gCreate a file to store hostname allocated
export HOSTFILE=$SLURM_SUBMIT_DIR/host-$SLURM_JOBID
module purge

module load firedrake/20240327

srun hostname > $HOSTFILE
## Information about the entry and exit of the job
echo -e "\n## Diretorio de submissao do job: $SLURM_SUBMIT_DIR \n"

# mpiexec -n 1 python save_new_shot_records.py
# mpiexec -n 12 python shot_filters.py
mpiexec -n 40 python test_small_marmousi.py 4 5.0 4.0

echo -e "\n## Job finished on $(date +'%d-%m-%Y as %T') ###################"
rm $HOSTFILE

rm -rf $FIREDRAKE_CACHE_DIR
25 changes: 24 additions & 1 deletion spyro/io/model_parameters.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,29 @@
from scipy.signal import sosfilt, iirfilter, zpk2sos


def cells_per_wavelength(method, degree, dimension):
if method == "mass_lumped_triangle":
method = "MLT"
cell_per_wavelength_dictionary = {
'mlt2tri': 7.02,
'mlt3tri': 3.70,
'mlt4tri': 2.67,
'mlt5tri': 2.03,
'mlt2tet': 6.12,
'mlt3tet': 3.72,
}
print(f"method: {method}", flush=True)

if dimension == 2 and (method == 'MLT' or method == 'CG'):
cell_type = 'tri'
if dimension == 3 and (method == 'MLT' or method == 'CG'):
cell_type = 'tet'

key = method.lower()+str(degree)+cell_type

return cell_per_wavelength_dictionary.get(key)


def butter_lowpass_filter_source(wavelet, cutoff, fs, order=2):
"""Low-pass filter the shot record with sampling-rate fs Hz
and cutoff freq. Hz
Expand Down Expand Up @@ -775,7 +798,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)
mesh_parameters.setdefault("cells_per_wavelength", None)
mesh_parameters.setdefault("cells_per_wavelength", cells_per_wavelength(self.method, self.degree, self.dimension))

self._set_mesh_length(
length_z=mesh_parameters["length_z"],
Expand Down
19 changes: 11 additions & 8 deletions spyro/meshing/meshing_functions.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,9 @@ def cells_per_wavelength(method, degree, dimension):
'mlt3tet': 3.72,
}

if dimension == 2 and (method == 'KMV' or method == 'CG'):
if dimension == 2 and (method == 'MLT' or method == 'CG'):
cell_type = 'tri'
if dimension == 3 and (method == 'KMV' or method == 'CG'):
if dimension == 3 and (method == 'MLT' or method == 'CG'):
cell_type = 'tet'

key = method.lower()+str(degree)+cell_type
Expand Down Expand Up @@ -358,6 +358,7 @@ def create_seimicmesh_2d_mesh(self):
"""
Creates a 2D mesh based on SeismicMesh meshing utilities.
"""
print(f"velocity_model{self.velocity_model}", flush=True)
if self.velocity_model is None:
return self.create_seismicmesh_2D_mesh_homogeneous()
else:
Expand All @@ -367,17 +368,19 @@ def create_seismicmesh_2D_mesh_with_velocity_model(self):
if self.comm.ensemble_comm.rank == 0:
v_min = self.minimum_velocity
frequency = self.source_frequency
if self.cpw is None:
self.cpw = cells_per_wavelength('MLT', 4, 2)
C = self.cpw # cells_per_wavelength(method, degree, dimension)

Lz = self.length_z*1000
Lx = self.length_x*1000
domain_pad = self.abc_pad*1000
Lz = self.length_z
Lx = self.length_x
domain_pad = self.abc_pad
lbda_min = v_min/frequency

bbox = (-Lz, 0.0, 0.0, Lx)
domain = SeismicMesh.Rectangle(bbox)

hmin = lbda_min/C*1000
hmin = lbda_min/C
self.comm.comm.barrier()

ef = SeismicMesh.get_sizing_function_from_segy(
Expand Down Expand Up @@ -408,15 +411,15 @@ def create_seismicmesh_2D_mesh_with_velocity_model(self):
if self.comm.comm.rank == 0:
meshio.write_points_cells(
"automatic_mesh.msh",
points/1000.0,
points,
[("triangle", cells)],
file_format="gmsh22",
binary=False
)

meshio.write_points_cells(
"automatic_mesh.vtk",
points/1000.0,
points,
[("triangle", cells)],
file_format="vtk"
)
Expand Down
1 change: 1 addition & 0 deletions spyro/solvers/inversion.py
Original file line number Diff line number Diff line change
Expand Up @@ -219,6 +219,7 @@ def generate_real_shot_record(self, plot_model=False, model_filename="model.png"
if Wave_obj_real_velocity.mesh is None and self.real_mesh is not None:
Wave_obj_real_velocity.mesh = self.real_mesh
if Wave_obj_real_velocity.initial_velocity_model is None:
print("C", flush=True)
Wave_obj_real_velocity.initial_velocity_model = self.real_velocity_model

if plot_model and Wave_obj_real_velocity.comm.comm.rank == 0 and Wave_obj_real_velocity.comm.ensemble_comm.rank == 0:
Expand Down
4 changes: 3 additions & 1 deletion spyro/solvers/wave.py
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ def matrix_building(self):
def set_mesh(
self,
user_mesh=None,
mesh_parameters=None,
mesh_parameters={},
):
"""
Set the mesh for the solver.
Expand Down Expand Up @@ -170,6 +170,8 @@ def set_initial_velocity_model(
output: bool (optional)
If True, outputs the velocity model to a pvd file for visualization.
"""
if new_file is not None:
self.initial_velocity_model_file = new_file
# If no mesh is set, we have to do it beforehand
if self.mesh is None:
self.set_mesh()
Expand Down
89 changes: 89 additions & 0 deletions test_small_marmousi.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,89 @@
import spyro
import sys


degree = int(sys.argv[1])
frequency = float(sys.argv[2])
final_time = float(sys.argv[3])


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": degree, # p order
"dimension": 2, # dimension
}
dictionary["parallelism"] = {
"type": "automatic", # options: automatic (same number of cores for evey processor) or spatial
}
dictionary["mesh"] = {
"Lz": 3.5, # depth in km - always positive # Como ver isso sem ler a malha?
"Lx": 17.0, # width in km - always positive
"Ly": 0.0, # thickness in km - always positive
"mesh_file": None,
}
dictionary["acquisition"] = {
"source_type": "ricker",
"source_locations": spyro.create_transect((-0.01, 4.0), (-0.01, 12.0), 20),
# "source_locations": [(-0.01, 4.0)],
"frequency": frequency,
# "frequency_filter": frequency_filter,
"delay": 0.2,
"delay_type": "time",
"receiver_locations": spyro.create_transect((-1.45, 4.0), (-1.45, 12.0), 100),
}
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
"gradient_sampling_frequency": 1, # how frequently to save solution to RAM
}
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():
dictionary["mesh"]["mesh_file"] = "meshes/real5hz.msh"

fwi = spyro.FullWaveformInversion(dictionary=dictionary)
fwi.set_real_velocity_model(new_file="velocity_models/vp_marmousi-ii.hdf5")
fwi.generate_real_shot_record(plot_model=True, save_shot_record=True)


def test_realistic_fwi():
dictionary["inversion"] = {
"perform_fwi": True,
"real_shot_record_files": f"shots/shot_record_",
"initial_guess_model_file": "velocity_models/initial_guess_15Hz.hdf5",
}
fwi = spyro.FullWaveformInversion(dictionary=dictionary)
fwi.set_guess_velocity_model(new_file="velocity_models/initial_guess_15Hz.hdf5")
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.5, maxiter=60)


if __name__ == "__main__":
test_real_shot_record_generation_parallel()
# test_realistic_fwi()

0 comments on commit 9c93b49

Please sign in to comment.