From f9458bab22772fb125ec0adfee1ec309e2d7703f Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 18 Mar 2024 07:32:43 -0400 Subject: [PATCH] remove duplicate so --- src/TimeModeling/Modeling/misfit_fg.jl | 2 +- src/TimeModeling/Modeling/python_interface.jl | 20 ++--- src/pysource/FD_utils.py | 13 +-- src/pysource/interface.py | 86 +++++++------------ src/pysource/propagators.py | 16 +++- src/pysource/sensitivity.py | 7 +- 6 files changed, 60 insertions(+), 84 deletions(-) diff --git a/src/TimeModeling/Modeling/misfit_fg.jl b/src/TimeModeling/Modeling/misfit_fg.jl index 4ad179cb7..585e21b01 100644 --- a/src/TimeModeling/Modeling/misfit_fg.jl +++ b/src/TimeModeling/Modeling/misfit_fg.jl @@ -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) diff --git a/src/TimeModeling/Modeling/python_interface.jl b/src/TimeModeling/Modeling/python_interface.jl index 50c42d636..7e49d9d3a 100644 --- a/src/TimeModeling/Modeling/python_interface.jl +++ b/src/TimeModeling/Modeling/python_interface.jl @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/pysource/FD_utils.py b/src/pysource/FD_utils.py index 81738952b..2dd7977cf 100644 --- a/src/pysource/FD_utils.py +++ b/src/pysource/FD_utils.py @@ -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 @@ -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) diff --git a/src/pysource/interface.py b/src/pysource/interface.py index 319a8b0eb..0b00945f4 100644 --- a/src/pysource/interface.py +++ b/src/pysource/interface.py @@ -14,7 +14,7 @@ # Forward wrappers Pr*F*Ps'*q -def forward_rec(model, src_coords, wavelet, rec_coords, space_order=8, f0=0.015, +def forward_rec(model, src_coords, wavelet, rec_coords, f0=0.015, illum=False, fw=True): """ Modeling of a point source with receivers Pr*F*Ps^T*q. @@ -29,8 +29,6 @@ def forward_rec(model, src_coords, wavelet, rec_coords, space_order=8, f0=0.015, Source signature rec_coords: Array Coordiantes of the receiver(s) - space_order: Int (optional) - Spatial discretization order, defaults to 8 f0: float peak frequency illum: bool @@ -44,12 +42,12 @@ def forward_rec(model, src_coords, wavelet, rec_coords, space_order=8, f0=0.015, Shot record """ rec, _, I, _ = forward(model, src_coords, rec_coords, wavelet, save=False, - space_order=space_order, f0=f0, illum=illum, fw=fw) + f0=f0, illum=illum, fw=fw) return rec.data, getattr(I, "data", None) # Pr*F*Pw'*w -def forward_rec_w(model, weight, wavelet, rec_coords, space_order=8, f0=0.015, +def forward_rec_w(model, weight, wavelet, rec_coords, f0=0.015, illum=False, fw=True): """ Forward modeling of an extended source with receivers Pr*F*Pw^T*w @@ -64,8 +62,6 @@ def forward_rec_w(model, weight, wavelet, rec_coords, space_order=8, f0=0.015, Source signature rec_coords: Array Coordiantes of the receiver(s) - space_order: Int (optional) - Spatial discretization order, defaults to 8 f0: float peak frequency illum: bool @@ -79,12 +75,12 @@ def forward_rec_w(model, weight, wavelet, rec_coords, space_order=8, f0=0.015, Shot record """ rec, _, I, _ = forward(model, None, rec_coords, wavelet, save=False, ws=weight, - space_order=space_order, f0=f0, illum=illum, fw=fw) + f0=f0, illum=illum, fw=fw) return rec.data, getattr(I, "data", None) # F*Ps'*q -def forward_no_rec(model, src_coords, wavelet, space_order=8, f0=0.015, illum=False, +def forward_no_rec(model, src_coords, wavelet, f0=0.015, illum=False, fw=True): """ Forward modeling of a point source without receiver. @@ -97,8 +93,6 @@ def forward_no_rec(model, src_coords, wavelet, space_order=8, f0=0.015, illum=Fa Coordiantes of the source(s) wavelet: Array Source signature - space_order: Int (optional) - Spatial discretization order, defaults to 8 f0: float peak frequency illum: bool @@ -111,13 +105,13 @@ def forward_no_rec(model, src_coords, wavelet, space_order=8, f0=0.015, illum=Fa Array Wavefield """ - _, u, I, _ = forward(model, src_coords, None, wavelet, space_order=space_order, + _, u, I, _ = forward(model, src_coords, None, wavelet, save=True, f0=f0, illum=illum, fw=fw) return u.data, getattr(I, "data", None) # Pr*F*u -def forward_wf_src(model, u, rec_coords, space_order=8, f0=0.015, illum=False, fw=True): +def forward_wf_src(model, u, rec_coords, f0=0.015, illum=False, fw=True): """ Forward modeling of a full wavefield source Pr*F*u. @@ -129,8 +123,6 @@ def forward_wf_src(model, u, rec_coords, space_order=8, f0=0.015, illum=False, f Time-space dependent wavefield rec_coords: Array Coordiantes of the receiver(s) - space_order: Int (optional) - Spatial discretization order, defaults to 8 f0: float peak frequency illum: bool @@ -144,13 +136,13 @@ def forward_wf_src(model, u, rec_coords, space_order=8, f0=0.015, illum=False, f Shot record """ wsrc = src_wavefield(model, u, fw=True) - rec, _, I, _ = forward(model, None, rec_coords, None, space_order=space_order, + rec, _, I, _ = forward(model, None, rec_coords, None, qwf=wsrc, illum=illum, f0=f0, fw=fw) return rec.data, getattr(I, "data", None) # F*u -def forward_wf_src_norec(model, u, space_order=8, f0=0.015, illum=False, fw=True): +def forward_wf_src_norec(model, u, f0=0.015, illum=False, fw=True): """ Forward modeling of a full wavefield source without receiver F*u. @@ -160,8 +152,6 @@ def forward_wf_src_norec(model, u, space_order=8, f0=0.015, illum=False, fw=True Physical model u: TimeFunction or Array Time-space dependent wavefield - space_order: Int (optional) - Spatial discretization order, defaults to 8 f0: float peak frequency illum: bool @@ -175,13 +165,13 @@ def forward_wf_src_norec(model, u, space_order=8, f0=0.015, illum=False, fw=True Wavefield """ wf_src = src_wavefield(model, u, fw=True) - _, u, I, _ = forward(model, None, None, None, space_order=space_order, save=True, + _, u, I, _ = forward(model, None, None, None, save=True, qwf=wf_src, f0=f0, illum=illum, fw=fw) return u.data, getattr(I, "data", None) # Pw*F'*Pr'*d_obs -def adjoint_w(model, rec_coords, data, wavelet, space_order=8, f0=0.015, illum=False, +def adjoint_w(model, rec_coords, data, wavelet, f0=0.015, illum=False, fw=True): """ Adjoint/backward modeling of a shot record (receivers as source) for an @@ -197,8 +187,6 @@ def adjoint_w(model, rec_coords, data, wavelet, space_order=8, f0=0.015, illum=F Shot gather wavelet: Array Time signature of the forward source for stacking along time - space_order: Int (optional) - Spatial discretization order, defaults to 8 f0: float peak frequency illum: bool @@ -212,13 +200,13 @@ def adjoint_w(model, rec_coords, data, wavelet, space_order=8, f0=0.015, illum=F spatial distribution """ w, _, I, _ = forward(model, rec_coords, None, data, wr=wavelet, - space_order=space_order, f0=f0, illum=illum, fw=fw) + f0=f0, illum=illum, fw=fw) return w.data, getattr(I, "data", None) # Linearized modeling ∂/∂m (Pr*F*Ps'*q) def born_rec(model, src_coords, wavelet, rec_coords, - space_order=8, ic="as", f0=0.015, illum=False, fw=True): + ic="as", f0=0.015, illum=False, fw=True): """ Linearized (Born) modeling of a point source for a model perturbation (square slowness) dm. @@ -233,8 +221,6 @@ def born_rec(model, src_coords, wavelet, rec_coords, Source signature rec_coords: Array Coordiantes of the receiver(s) - space_order: Int (optional) - Spatial discretization order, defaults to 8 ic: String Imaging conditions ("as", "isic" or "fwi"), defaults to "as" f0: float @@ -250,13 +236,13 @@ def born_rec(model, src_coords, wavelet, rec_coords, Shot record """ rec, _, I, _ = born(model, src_coords, rec_coords, wavelet, save=False, - space_order=space_order, ic=ic, f0=f0, illum=illum, fw=fw) + ic=ic, f0=f0, illum=illum, fw=fw) return rec.data, getattr(I, "data", None) # ∂/∂m (Pr*F*Pw'*w) def born_rec_w(model, weight, wavelet, rec_coords, - space_order=8, ic="as", f0=0.015, illum=False, fw=True): + ic="as", f0=0.015, illum=False, fw=True): """ Linearized (Born) modeling of an extended source for a model perturbation (square slowness) dm with an extended source @@ -271,8 +257,6 @@ def born_rec_w(model, weight, wavelet, rec_coords, Source signature rec_coords: Array Coordiantes of the receiver(s) - space_order: Int (optional) - Spatial discretization order, defaults to 8 ic: String Imaging conditions ("as", "isic" or "fwi"), defaults to "as" f0: float @@ -288,11 +272,11 @@ def born_rec_w(model, weight, wavelet, rec_coords, Shot record """ rec, _, I, _ = born(model, None, rec_coords, wavelet, save=False, ws=weight, - space_order=space_order, ic=ic, f0=f0, illum=illum, fw=fw) + ic=ic, f0=f0, illum=illum, fw=fw) return rec.data, getattr(I, "data", None) -def J_adjoint(model, src_coords, wavelet, rec_coords, recin, space_order=8, +def J_adjoint(model, src_coords, wavelet, rec_coords, recin, is_residual=False, checkpointing=False, n_checkpoints=None, t_sub=1, return_obj=False, freq_list=[], dft_sub=None, ic="as", illum=False, ws=None, f0=0.015, born_fwd=False, nlind=False, misfit=None, fw=True): @@ -315,8 +299,6 @@ def J_adjoint(model, src_coords, wavelet, rec_coords, recin, space_order=8, Coordiantes of the receiver(s) recin: Array Receiver data - space_order: Int (optional) - Spatial discretization order, defaults to 8 checkpointing: Bool Whether or not to use checkpointing n_checkpoints: Int @@ -345,25 +327,25 @@ def J_adjoint(model, src_coords, wavelet, rec_coords, recin, space_order=8, """ if checkpointing: return J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, - space_order=8, is_residual=is_residual, ws=ws, + is_residual=is_residual, ws=ws, n_checkpoints=n_checkpoints, ic=ic, f0=f0, nlind=nlind, return_obj=return_obj, illum=illum, born_fwd=born_fwd, misfit=misfit, fw=fw) elif freq_list is not None: return J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, ws=ws, - space_order=space_order, dft_sub=dft_sub, f0=f0, ic=ic, + dft_sub=dft_sub, f0=f0, ic=ic, freq_list=freq_list, is_residual=is_residual, nlind=nlind, return_obj=return_obj, misfit=misfit, born_fwd=born_fwd, illum=illum, fw=fw) else: return J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, is_residual=is_residual, ic=ic, ws=ws, t_sub=t_sub, - return_obj=return_obj, space_order=space_order, + return_obj=return_obj, born_fwd=born_fwd, f0=f0, nlind=nlind, illum=illum, misfit=misfit, fw=fw) -def J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, space_order=8, +def J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, freq_list=[], is_residual=False, return_obj=False, nlind=False, dft_sub=None, ic="as", ws=None, born_fwd=False, f0=0.015, misfit=None, illum=False, fw=True): @@ -384,8 +366,6 @@ def J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, space_order=8, Coordiantes of the receiver(s) recin: Array Receiver data - space_order: Int (optional) - Spatial discretization order, defaults to 8 freq_list: List List of frequencies for on-the-fly DFT dft_sub: Int @@ -415,20 +395,20 @@ def J_adjoint_freq(model, src_coords, wavelet, rec_coords, recin, space_order=8, """ ffunc = op_fwd_J[born_fwd] rec, u, Iu, _ = ffunc(model, src_coords, rec_coords, wavelet, save=False, - space_order=space_order, freq_list=freq_list, ic=ic, ws=ws, + freq_list=freq_list, ic=ic, ws=ws, dft_sub=dft_sub, nlind=nlind, illum=illum, f0=f0, fw=fw) # Residual and gradient f, residual = Loss(rec, recin, model.critical_dt, is_residual=is_residual, misfit=misfit) - g, Iv, _ = gradient(model, residual, rec_coords, u, space_order=space_order, ic=ic, + g, Iv, _ = gradient(model, residual, rec_coords, u, ic=ic, freq=freq_list, dft_sub=dft_sub, f0=f0, illum=illum, fw=fw) if return_obj: return f, g.data, getattr(Iu, "data", None), getattr(Iv, "data", None) return g.data, getattr(Iu, "data", None), getattr(Iv, "data", None) -def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, space_order=8, +def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, is_residual=False, return_obj=False, born_fwd=False, illum=False, ic="as", ws=None, t_sub=1, nlind=False, f0=0.015, misfit=None, fw=True): @@ -449,8 +429,6 @@ def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, space_orde Coordiantes of the receiver(s) recin: Array Receiver data - space_order: Int (optional) - Spatial discretization order, defaults to 8 ic: String Imaging conditions ("as", "isic" or "fwi"), defaults to "as" ws : Array @@ -476,14 +454,14 @@ def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, space_orde """ ffunc = op_fwd_J[born_fwd] rec, u, Iu, _ = ffunc(model, src_coords, rec_coords, wavelet, save=True, nlind=nlind, - f0=f0, ws=ws, space_order=space_order, illum=illum, ic=ic, + f0=f0, ws=ws, illum=illum, ic=ic, t_sub=t_sub, fw=fw) # Residual and gradient f, residual = Loss(rec, recin, model.critical_dt, is_residual=is_residual, misfit=misfit) - g, Iv, _ = gradient(model, residual, rec_coords, u, space_order=space_order, ic=ic, + g, Iv, _ = gradient(model, residual, rec_coords, u, ic=ic, f0=f0, illum=illum, fw=fw) if return_obj: @@ -491,7 +469,7 @@ def J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin, space_orde return g.data, getattr(Iu, "data", None), getattr(Iv, "data", None) -def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space_order=8, +def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, is_residual=False, n_checkpoints=None, born_fwd=False, return_obj=False, ic="as", ws=None, nlind=False, f0=0.015, misfit=None, illum=False, fw=True): @@ -511,8 +489,6 @@ def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space Coordiantes of the receiver(s) recin: Array Receiver data - space_order: Int (optional) - Spatial discretization order, defaults to 8 checkpointing: Bool Whether or not to use checkpointing n_checkpoints: Int @@ -545,9 +521,9 @@ def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space ffunc = op_fwd_J[born_fwd] # Optimal checkpointing op_f, u, rec_g, kwu = ffunc(model, src_coords, rec_coords, wavelet, fw=fw, - save=False, space_order=space_order, return_op=True, + save=False, return_op=True, ic=ic, nlind=nlind, ws=ws, f0=f0, illum=illum) - op, g, kwg = gradient(model, recin, rec_coords, u, space_order=space_order, + op, g, kwg = gradient(model, recin, rec_coords, u, return_op=True, ic=ic, f0=f0, save=False, illum=illum, fw=fw) @@ -588,7 +564,7 @@ def J_adjoint_checkpointing(model, src_coords, wavelet, rec_coords, recin, space op_fwd_J = {False: forward, True: born} -def wri_func(model, src_coords, wavelet, rec_coords, recin, yin, space_order=8, +def wri_func(model, src_coords, wavelet, rec_coords, recin, yin, ic="as", ws=None, t_sub=1, grad="m", grad_corr=False, alpha_op=False, w_fun=None, eps=0, freq_list=[], wfilt=None, f0=0.015): """ @@ -606,7 +582,7 @@ def wri_func(model, src_coords, wavelet, rec_coords, recin, yin, space_order=8, # F(m0) * q if y is not an input and compute y = r(m0) if yin is None or grad_corr: y, u0, _, _ = forward(model, src_coords, rec_coords, wavelet, save=grad_corr, - space_order=space_order, ws=ws, f0=f0) + ws=ws, f0=f0) ydat = recin[:] - y.data[:] else: ydat = yin diff --git a/src/pysource/propagators.py b/src/pysource/propagators.py index 4481df0bd..dbaf59681 100644 --- a/src/pysource/propagators.py +++ b/src/pysource/propagators.py @@ -12,7 +12,7 @@ # Forward propagation -def forward(model, src_coords, rcv_coords, wavelet, space_order=8, save=False, +def forward(model, src_coords, rcv_coords, wavelet, save=False, qwf=None, return_op=False, freq_list=None, dft_sub=None, norm_wf=False, w_fun=None, ws=None, wr=None, t_sub=1, f0=0.015, illum=False, fw=True, **kwargs): @@ -20,6 +20,8 @@ def forward(model, src_coords, rcv_coords, wavelet, space_order=8, save=False, Low level propagator, to be used through `interface.py` Compute forward wavefield u = A(m)^{-1}*f and related quantities (u(xrcv)) """ + # Space order + space_order = model.space_order # Number of time steps nt = as_tuple(qwf)[0].shape[0] if wavelet is None else wavelet.shape[0] @@ -92,12 +94,14 @@ def adjoint(*args, **kwargs): return forward(*args, fw=fw, **kwargs) -def gradient(model, residual, rcv_coords, u, return_op=False, space_order=8, fw=True, +def gradient(model, residual, rcv_coords, u, return_op=False, fw=True, w=None, freq=None, dft_sub=None, ic="as", f0=0.015, save=True, illum=False): """ Low level propagator, to be used through `interface.py` Compute the action of the adjoint Jacobian onto a residual J'* δ d. """ + # Space order + space_order = model.space_order # Setting adjoint wavefieldgradient v = wavefield(model, space_order, fw=not fw) try: @@ -144,7 +148,7 @@ def gradient(model, residual, rcv_coords, u, return_op=False, space_order=8, fw= return gradm, I, summary -def born(model, src_coords, rcv_coords, wavelet, space_order=8, save=False, +def born(model, src_coords, rcv_coords, wavelet, save=False, qwf=None, return_op=False, ic="as", freq_list=None, dft_sub=None, ws=None, t_sub=1, nlind=False, f0=0.015, illum=False, fw=True): """ @@ -152,6 +156,8 @@ def born(model, src_coords, rcv_coords, wavelet, space_order=8, save=False, Compute linearized wavefield U = J(m)* δ m and related quantities. """ + # Space order + space_order = model.space_order nt = wavelet.shape[0] # Wavefields @@ -219,12 +225,14 @@ def born(model, src_coords, rcv_coords, wavelet, space_order=8, save=False, # Forward propagation -def forward_grad(model, src_coords, rcv_coords, wavelet, v, space_order=8, +def forward_grad(model, src_coords, rcv_coords, wavelet, v, q=None, ws=None, ic="as", w=None, freq=None, f0=0.015, **kwargs): """ Low level propagator, to be used through `interface.py` Compute forward wavefield u = A(m)^{-1}*f and related quantities (u(xrcv)) """ + # Space order + space_order = model.space_order # Number of time steps nt = as_tuple(q)[0].shape[0] if wavelet is None else wavelet.shape[0] diff --git a/src/pysource/sensitivity.py b/src/pysource/sensitivity.py index fbeb3f903..e06a0fbd8 100644 --- a/src/pysource/sensitivity.py +++ b/src/pysource/sensitivity.py @@ -205,8 +205,7 @@ def isic_src(model, u, **kwargs): ics = kwargs.get('icsign', 1) dus = [] for uu in u: - so = uu.space_order // 2 - du_aux = div(dm * irho * grad(uu, shift=.5, order=so), shift=-.5, order=so) + du_aux = div(dm * irho * grad(uu, shift=.5), shift=-.5) dus.append(dm * irho * uu.dt2 * m - ics * du_aux) if model.is_tti: return (-dus[0], -dus[1]) @@ -224,9 +223,7 @@ def inner_grad(u, v): v: TimeFunction or Function Second wavefield """ - so = u.space_order // 2 - return sum([a*b for a, b in zip(grad(u, order=so), - grad(v, order=so))]) + return sum([a*b for a, b in zip(grad(u), grad(v))]) fwi_src = lambda *ar, **kw: isic_src(*ar, icsign=-1, **kw)