Skip to content

Commit

Permalink
remove duplicate so
Browse files Browse the repository at this point in the history
  • Loading branch information
mloubout committed Mar 18, 2024
1 parent 7fd4d13 commit f9458ba
Show file tree
Hide file tree
Showing 6 changed files with 60 additions and 84 deletions.
2 changes: 1 addition & 1 deletion src/TimeModeling/Modeling/misfit_fg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ function multi_src_fg(model_full::AbstractModel, source::judiVector, dObs::judiV
@juditime "Python call to J_adjoint" begin
argout = rlock_pycall(ac."J_adjoint", Tuple{Float32, PyArray, IT...}, modelPy,
src_coords, qIn, rec_coords, dObserved, t_sub=options.subsampling_factor,
space_order=options.space_order, checkpointing=options.optimal_checkpointing,
checkpointing=options.optimal_checkpointing,
freq_list=freqs, ic=options.IC, is_residual=false, born_fwd=lin, nlind=nlind,
dft_sub=options.dft_subsampling_factor[1], f0=options.f0, return_obj=true,
misfit=mfunc, illum=illum)
Expand Down
20 changes: 10 additions & 10 deletions src/TimeModeling/Modeling/python_interface.jl
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Geometry, srcData::Arr
rec_coords = setup_grid(recGeometry, modelPy.shape)

# Devito call
return wrapcall_data(ac."forward_rec", modelPy, src_coords, qIn, rec_coords, fw=fw, space_order=options.space_order, f0=options.f0, illum=illum)
return wrapcall_data(ac."forward_rec", modelPy, src_coords, qIn, rec_coords, fw=fw, f0=options.f0, illum=illum)
end

# u = F*Ps'*q
Expand All @@ -76,7 +76,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Geometry, srcData::Arr
src_coords = setup_grid(srcGeometry, modelPy.shape)

# Devito call
return wrapcall_wf(ac."forward_no_rec", modelPy, src_coords, qIn, fw=fw, space_order=options.space_order, illum=illum)
return wrapcall_wf(ac."forward_no_rec", modelPy, src_coords, qIn, fw=fw, illum=illum)
end

# d_obs = Pr*F*u
Expand All @@ -89,7 +89,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Nothing, srcData::Arra
rec_coords = setup_grid(recGeometry, modelPy.shape)

# Devito call
return wrapcall_data(ac."forward_wf_src", modelPy, srcData, rec_coords, fw=fw, space_order=options.space_order, f0=options.f0, illum=illum)
return wrapcall_data(ac."forward_wf_src", modelPy, srcData, rec_coords, fw=fw, f0=options.f0, illum=illum)
end

# u_out = F*u_in
Expand All @@ -99,7 +99,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Nothing, srcData::Arra
dtComp = convert(Float32, modelPy."critical_dt")

# Devito call
return wrapcall_wf(ac."forward_wf_src_norec", modelPy, srcData, fw=fw, space_order=options.space_order, illum=illum)
return wrapcall_wf(ac."forward_wf_src_norec", modelPy, srcData, fw=fw, illum=illum)
end

# d_lin = J*dm
Expand All @@ -117,7 +117,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Geometry, srcData::Arr

# Devito call
return wrapcall_data(ac."born_rec", modelPy, src_coords, qIn, rec_coords, fw=fw,
space_order=options.space_order, ic=options.IC, f0=options.f0, illum=illum)
ic=options.IC, f0=options.f0, illum=illum)
end

# dm = J'*d_lin
Expand All @@ -136,7 +136,7 @@ function devito_interface(modelPy::PyObject, srcGeometry::Geometry, srcData::Arr
length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies
return wrapcall_grad(ac."J_adjoint", modelPy,
src_coords, qIn, rec_coords, dIn, fw=fw, t_sub=options.subsampling_factor,
space_order=options.space_order, checkpointing=options.optimal_checkpointing,
checkpointing=options.optimal_checkpointing,
freq_list=freqs, ic=options.IC, is_residual=true,
dft_sub=options.dft_subsampling_factor[1], f0=options.f0, illum=illum)
end
Expand All @@ -157,7 +157,7 @@ function devito_interface(modelPy::PyObject, weights::Array, srcData::Array, rec

# Devito call
return wrapcall_data(ac."forward_rec_w", modelPy, weights, qIn, rec_coords,
fw=fw, space_order=options.space_order, f0=options.f0, illum=illum)
fw=fw, f0=options.f0, illum=illum)
end

# dw = Pw*F'*Pr'*d_obs - adjoint modeling w/ extended source
Expand All @@ -173,7 +173,7 @@ function devito_interface(modelPy::PyObject, recGeometry::Geometry, recData::Arr

# Devito call
return wrapcall_weights(ac."adjoint_w", modelPy, rec_coords, dIn, qIn,
fw=fw, space_order=options.space_order, f0=options.f0, illum=illum)
fw=fw, f0=options.f0, illum=illum)
end

# Jacobian of extended source modeling: d_lin = J*dm
Expand All @@ -190,7 +190,7 @@ function devito_interface(modelPy::PyObject, weights::Array, srcData::Array, rec

# Devito call
return wrapcall_data(ac."born_rec_w", modelPy, weights, qIn, rec_coords,
fw=fw, space_order=options.space_order, ic=options.IC, f0=options.f0, illum=illum)
fw=fw, ic=options.IC, f0=options.f0, illum=illum)
end

# Adjoint Jacobian of extended source modeling: dm = J'*d_lin
Expand All @@ -208,7 +208,7 @@ function devito_interface(modelPy::PyObject, weights::Array, srcData::Array, rec
length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies
return wrapcall_grad(ac."J_adjoint", modelPy,
nothing, qIn, rec_coords, dIn, fw=fw, t_sub=options.subsampling_factor,
space_order=options.space_order, checkpointing=options.optimal_checkpointing,
checkpointing=options.optimal_checkpointing,
freq_list=freqs, ic=options.IC, ws=weights, is_residual=true,
dft_sub=options.dft_subsampling_factor[1], f0=options.f0, illum=illum)
end
13 changes: 4 additions & 9 deletions src/pysource/FD_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,7 @@ def laplacian(v, irho):
"""
irho = irho or 1
if isinstance(irho, Differentiable):
so = irho.space_order // 2
Lap = div(irho * grad(v, shift=.5, order=so), shift=-.5, order=so)
Lap = div(irho * grad(v, shift=.5), shift=-.5)
else:
Lap = irho * v.laplace

Expand Down Expand Up @@ -93,16 +92,12 @@ def sa_tti(u, v, model):
model: Model
Model structure
"""
# Space order
so = u.space_order // 2
# Matrix of Thomsen params
A, B, C = thomsen_mat(model)
# Rotation Matrix
R = R_mat(model)

PI = R.T * (A * R * grad(u, shift=.5, order=so) +
B * R * grad(v, shift=.5, order=so))
MI = R.T * (B * R * grad(u, shift=.5, order=so) +
C * R * grad(v, shift=.5, order=so))
PI = R.T * (A * R * grad(u, shift=.5) + B * R * grad(v, shift=.5))
MI = R.T * (B * R * grad(u, shift=.5) + C * R * grad(v, shift=.5))

return div(PI, shift=-.5, order=so), div(MI, shift=-.5, order=so)
return div(PI, shift=-.5), div(MI, shift=-.5)
Loading

0 comments on commit f9458ba

Please sign in to comment.