Skip to content

Commit

Permalink
small changes just to update dictionary to new model
Browse files Browse the repository at this point in the history
  • Loading branch information
Olender committed Nov 27, 2024
1 parent f013d58 commit 1651eb7
Show file tree
Hide file tree
Showing 3 changed files with 26 additions and 47 deletions.
17 changes: 4 additions & 13 deletions spyro/solvers/forward_ad.py
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ def execute_acoustic(
# Sources.
source_mesh = fire.VertexOnlyMesh(
self.mesh,
[self.model["acquisition"]["source_pos"][source_number]]
[self.model["acquisition"]["source_locations"][source_number]]
)
# Source function space.
V_s = fire.FunctionSpace(source_mesh, "DG", 0)
Expand All @@ -88,7 +88,7 @@ def execute_acoustic(
# Time execution.
J_val = 0.0
receiver_data = []
total_steps = int(self.model["timeaxis"]["tf"] / self.model["timeaxis"]["dt"])
total_steps = int(self.model["time_axis"]["final_time"] / self.model["time_axis"]["dt"]) + 1
if (
fire_ad.get_working_tape()._checkpoint_manager
and self.model["aut_dif"]["checkpointing"]
Expand Down Expand Up @@ -119,19 +119,10 @@ def execute_elastic(self):
"for the automatic differentiation based FWI.")

def _solver_parameters(self):
if self.model["opts"]["method"] == "KMV":
if self.model["options"]["variant"] == "lumped":
params = {"ksp_type": "preonly", "pc_type": "jacobi"}
elif (
self.model["opts"]["method"] == "CG"
and self.mesh.ufl_cell() != quadrilateral # noqa: F821
and self.mesh.ufl_cell() != hexahedron # noqa: F821
):
elif self.model["options"]["variant"] == "equispaced":
params = {"ksp_type": "cg", "pc_type": "jacobi"}
elif self.model["opts"]["method"] == "CG" and (
self.mesh.ufl_cell() == quadrilateral # noqa: F821
or self.mesh.ufl_cell() == hexahedron # noqa: F821
):
params = {"ksp_type": "preonly", "pc_type": "jacobi"}
else:
raise ValueError("method is not yet supported")

Expand Down
4 changes: 2 additions & 2 deletions spyro/solvers/time_integration_ad.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ def central_difference_acoustic(forwardsolver, c, source_function):
"""
# Acoustic linear variational solver.
V = forwardsolver.V
dt = forwardsolver.model["timeaxis"]["dt"]
dt = forwardsolver.model["time_axis"]["dt"]
u = fire.TrialFunction(V)
v = fire.TestFunction(V)
u_np1 = fire.Function(V) # timestep n+1
Expand All @@ -30,7 +30,7 @@ def central_difference_acoustic(forwardsolver, c, source_function):
fire.Constant(dt**2) * v * fire.dx(scheme=qr_x)

nf = 0
if forwardsolver.model["BCs"]["outer_bc"] == "non-reflective":
if forwardsolver.model["absorving_boundary_conditions"]["status"] is True:
nf = (1/c) * ((u_n - u_nm1) / dt) * v * fire.ds(scheme=qr_s)

a = fire.dot(fire.grad(u_n), fire.grad(v)) * fire.dx(scheme=qr_x)
Expand Down
52 changes: 20 additions & 32 deletions test/test_gradient_ad.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,45 +8,33 @@
# --- Basid setup to run a forward simulation with AD --- #
model = {}

model["opts"] = {
"method": "KMV", # either CG or mass_lumped_triangle
"quadrature": "KMV", # Equi or mass_lumped_triangle
model["options"] = {
"cell_type": "T", # T or Q
"variant": "lumped",
"degree": 1, # p order
"dimension": 2, # dimension
"regularization": False, # regularization is on?
"gamma": 1e-5, # regularization parameter
}

model["parallelism"] = {
# options:
# `shots_parallelism`. Shots parallelism.
# None - no shots parallelism.
"type": None,
"num_spacial_cores": 1, # Number of cores to use in the spatial parallelism.
"type": "automatic", # options: automatic (same number of cores for evey processor) or spatial
}

# Define the domain size without the ABL.
model["mesh"] = {
"Lz": 1.0, # depth in km - always positive
"Lx": 1.0, # width in km - always positive
"Ly": 0.0, # thickness in km - always positive
"meshfile": "not_used.msh",
"initmodel": "not_used.hdf5",
"truemodel": "not_used.hdf5",
"mesh_file": None,
}

# Specify a 250-m Absorbing Boundary Layer (ABL) on the three sides of the domain to damp outgoing waves.
model["BCs"] = {
"status": False, # True or False, used to turn on any type of BC
"outer_bc": "non-reflective", # none or non-reflective (outer boundary condition)
"lz": 0.0, # thickness of the ABL in the z-direction (km) - always positive
"lx": 0.0, # thickness of the ABL in the x-direction (km) - always positive
"ly": 0.0, # thickness of the ABL in the y-direction (km) - always positive
model["absorving_boundary_conditions"] = {
"status": False,
"pad_length": 0.,
}

model["acquisition"] = {
"source_type": "Ricker",
"source_pos": spyro.create_transect((0.2, 0.15), (0.8, 0.15), 1),
"source_type": "ricker",
"source_locations": spyro.create_transect((0.2, 0.15), (0.8, 0.15), 1),
"frequency": 7.0,
"delay": 1.0,
"receiver_locations": spyro.create_transect((0.2, 0.2), (0.8, 0.2), 10),
Expand All @@ -56,13 +44,13 @@
"checkpointing": False,
}

model["timeaxis"] = {
"t0": 0.0, # Initial time for event
"tf": 0.4, # Final time for event (for test 7)
model["time_axis"] = {
"initial_time": 0.0, # Initial time for event
"final_time": 0.4, # Final time for event (for test 7)
"dt": 0.004, # timestep size (divided by 2 in the test 4. dt for test 3 is 0.00050)
"amplitude": 1, # the Ricker has an amplitude of 1.
"nspool": 20, # (20 for dt=0.00050) how frequently to output solution to pvds
"fspool": 1, # how frequently to save solution to RAM
"output_frequency": 20, # (20 for dt=0.00050) how frequently to output solution to pvds
"gradient_sampling_frequency": 1, # how frequently to save solution to RAM
}


Expand All @@ -89,7 +77,7 @@ def forward(
if annotate:
fire_ad.continue_annotation()
if model["aut_dif"]["checkpointing"]:
total_steps = int(model["timeaxis"]["tf"] / model["timeaxis"]["dt"])
total_steps = int(model["time_axis"]["final_time"] / model["time_axis"]["dt"]) + 1
steps_store = int(total_steps / 10) # Store 10% of the steps.
tape = fire_ad.get_working_tape()
tape.progress_bar = fire.ProgressBar
Expand All @@ -106,19 +94,19 @@ def forward(

def test_taylor():
# Test only serial for now.
M = model["parallelism"]["num_spacial_cores"]
M = 1
my_ensemble = fire.Ensemble(fire.COMM_WORLD, M)
mesh = fire.UnitSquareMesh(20, 20, comm=my_ensemble.comm)
element = fire.FiniteElement(
model["opts"]["method"], mesh.ufl_cell(), degree=model["opts"]["degree"],
variant=model["opts"]["quadrature"]
"KMV", mesh.ufl_cell(), degree=model["options"]["degree"],
variant="KMV"
)
V = fire.FunctionSpace(mesh, element)

fwd_solver = spyro.solvers.forward_ad.ForwardSolver(model, mesh, V)
# Ricker wavelet
wavelet = spyro.full_ricker_wavelet(
model["timeaxis"]["dt"], model["timeaxis"]["tf"],
model["time_axis"]["dt"], model["time_axis"]["final_time"],
model["acquisition"]["frequency"],
)
c_true = make_c_camembert(V, mesh)
Expand Down

0 comments on commit 1651eb7

Please sign in to comment.