Skip to content

Commit

Permalink
Merge pull request #95 from nogueirapeterson/viscoacoustic-operators
Browse files Browse the repository at this point in the history 
Including viscoacoustic operators
  • Loading branch information
@mloubout
mloubout authored Apr 20, 2022
2 parents cada5b3 + 1ed1f21 commit e2307a4
Show file tree
Hide file tree
Showing 31 changed files with 546 additions and 344 deletions.
51 changes: 33 additions & 18 deletions .github/workflows/ci-examples.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,41 +11,57 @@ on:
- master

jobs:
test:
name: JUDI example on Julia ${{ matrix.version }}
runs-on: ${{ matrix.os }}
list-examples:
runs-on: ubuntu-latest

strategy:
fail-fast: false

outputs:
matrix: ${{ steps.set-matrix.outputs.matrix }}

steps:
- name: Checkout JUDI
uses: actions/checkout@v2

- id: set-matrix
run: echo "::set-output name=matrix::$(ls examples/scripts/*.jl | xargs -n 1 basename | jq -R -s -c 'split("\n")[:-1]')"
shell: bash

run-examples:
runs-on: ubuntu-latest
needs: list-examples
name: ${{ matrix.example }} on Julia ${{ matrix.version }}

env:
DEVITO_ARCH: gcc-9
DEVITO_LANGUAGE: "openmp"
DEVITO_BACKEND: "core"
DEVITO_LOGGING: "ERROR"
OMP_NUM_THREADS: ${{ matrix.omp }}
OMP_NUM_THREADS: 2
NITER: 1

strategy:
fail-fast: false

matrix:
os: [ubuntu-latest]
omp: [2]
version: [1.6, 1.7]
example: ${{ fromJson(needs.list-examples.outputs.matrix) }}
version: ['1.7']

include:
- example: "modeling_basic_2D.jl"
version: '1.6'

steps:
- name: Checkout JUDI
uses: actions/checkout@v2

- name: Install GCC 9
if : runner.os == 'macOS'
run : brew install gcc@9

- name: Set julia python
run: |
PYTHON=$(which python3) julia -e 'using Pkg;Pkg.add("PyCall");Pkg.build("PyCall")'
- name: Setup julia
uses: julia-actions/setup-julia@v1
with:
version: "${{ matrix.version }}"
version: ${{ matrix.version }}
arch: x64

- name: Build JUDI
Expand All @@ -56,7 +72,6 @@ jobs:
julia -e 'using Pkg;Pkg.add(["NLopt", "JOLI", "PyPlot", "IterativeSolvers", "SlimOptim", "HDF5", "SegyIO", "SetIntersectionProjection"])'
julia -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd()))'
- name: Run base examples
working-directory: examples/scripts
run: |
julia -p 2 -e 'for f in filter!(e->e≠"runall.jl", readdir(".")); include(f); end'
- name: Run examples
working-directory: examples/scripts/
run: julia -p 2 ${{ matrix.example }}
12 changes: 11 additions & 1 deletion .github/workflows/ci-op.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ jobs:
matrix:
os: [ubuntu-latest]
op: ["ISO_OP", "ISO_OP_FS", "TTI_OP", "TTI_OP_FS", "BASICS"]
version: ['1.6', '1.7']
version: ['1.7']
omp: [2]

include:
Expand All @@ -41,6 +41,16 @@ jobs:
version: '1.7'
op: "ISO_OP_FS"
omp: 1

- os: ubuntu-latest
version: '1.7'
op: "VISCO_AC_OP"
omp: 2

- os: ubuntu-latest
version: '1.6'
op: "TTI_OP"
omp: 2

steps:
- name: Checkout JUDI
Expand Down
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "JUDI"
uuid = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
authors = ["Philipp Witte, Mathias Louboutin"]
version = "2.6.8"
version = "2.6.9"

This comment has been minimized.

Sign in to view

Copy link
@mloubout

mloubout Apr 20, 2022

Author Member

@JuliaRegistrator register

[deps]
DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
Expand Down
24 changes: 6 additions & 18 deletions src/TimeModeling/Modeling/fwi_objective_serial.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,25 +45,13 @@ function fwi_objective(model_full::Model, source::judiVector, dObs::judiVector,
src_coords = setup_grid(source.geometry, model.n) # shifts source coordinates by origin
rec_coords = setup_grid(dObs.geometry, model.n) # shifts rec coordinates by origin

length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies
argout1, argout2 = pycall(ac."J_adjoint", Tuple{Float32, PyArray}, modelPy,
src_coords, qIn, rec_coords, dObserved, t_sub=options.subsampling_factor,
space_order=options.space_order, checkpointing=options.optimal_checkpointing,
freq_list=freqs, isic=options.isic, is_residual=false, return_obj=true,
dft_sub=options.dft_subsampling_factor[1], f0=options.f0)

if options.optimal_checkpointing == true
argout1, argout2 = pycall(ac."J_adjoint_checkpointing", Tuple{Float32, PyArray},
modelPy, src_coords, qIn,
rec_coords, dObserved, is_residual=false, return_obj=true, isic=options.isic,
t_sub=options.subsampling_factor, space_order=options.space_order)
elseif ~isempty(options.frequencies)
argout1, argout2 = pycall(ac."J_adjoint_freq", Tuple{Float32, PyArray},
modelPy, src_coords, qIn,
rec_coords, dObserved, is_residual=false, return_obj=true, isic=options.isic,
freq_list=options.frequencies, t_sub=options.subsampling_factor,
space_order=options.space_order)
else
argout1, argout2 = pycall(ac."J_adjoint_standard", Tuple{Float32, PyArray},
modelPy, src_coords, qIn,
rec_coords, dObserved, is_residual=false, return_obj=true,
t_sub=options.subsampling_factor, space_order=options.space_order,
isic=options.isic)
end
argout2 = remove_padding(argout2, modelPy.padsizes; true_adjoint=options.sum_padding)
if options.limit_m==true
argout2 = extend_gradient(model_full, model, argout2)
Expand Down
26 changes: 7 additions & 19 deletions src/TimeModeling/Modeling/lsrtm_objective_serial.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,25 +48,13 @@ function lsrtm_objective(model_full::Model, source::judiVector, dObs::judiVector
src_coords = setup_grid(source.geometry, model.n) # shifts source coordinates by origin
rec_coords = setup_grid(dObs.geometry, model.n) # shifts rec coordinates by origin

if options.optimal_checkpointing == true
argout1, argout2 = pycall(ac."J_adjoint_checkpointing", Tuple{Float32, PyArray},
modelPy, src_coords, qIn,
rec_coords, dObserved, is_residual=false, return_obj=true,
t_sub=options.subsampling_factor, space_order=options.space_order,
born_fwd=true, nlind=nlind, isic=options.isic)
elseif ~isempty(options.frequencies)
argout1, argout2 = pycall(ac."J_adjoint_freq", Tuple{Float32, PyArray},
modelPy, src_coords, qIn,
rec_coords, dObserved, is_residual=false, return_obj=true, nlind=nlind,
freq_list=options.frequencies, t_sub=options.subsampling_factor,
space_order=options.space_order, born_fwd=true, isic=options.isic)
else
argout1, argout2 = pycall(ac."J_adjoint_standard", Tuple{Float32, PyArray},
modelPy, src_coords, qIn,
rec_coords, dObserved, is_residual=false, return_obj=true,
t_sub=options.subsampling_factor, space_order=options.space_order,
isic=options.isic, born_fwd=true, nlind=nlind)
end
length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies
argout1, argout2 = pycall(ac."J_adjoint", Tuple{Float32, PyArray}, modelPy,
src_coords, qIn, rec_coords, dObserved, t_sub=options.subsampling_factor,
space_order=options.space_order, checkpointing=options.optimal_checkpointing,
freq_list=freqs, isic=options.isic, is_residual=false, return_obj=true,
dft_sub=options.dft_subsampling_factor[1], f0=options.f0, born_fwd=true)

argout2 = remove_padding(argout2, modelPy.padsizes; true_adjoint=options.sum_padding)
if options.limit_m==true
argout2 = extend_gradient(model_full, model, argout2)
Expand Down
37 changes: 19 additions & 18 deletions src/TimeModeling/Modeling/python_interface.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Geometry, srcDa
rec_coords = setup_grid(recGeometry, modelPy.shape)

# Devito call
dOut = wrapcall_data(ac."forward_rec", modelPy, src_coords, qIn, rec_coords, space_order=options.space_order)
dOut = wrapcall_data(ac."forward_rec", modelPy, src_coords, qIn, rec_coords, space_order=options.space_order, f0=options.f0)
dOut = time_resample(dOut, dtComp, recGeometry)

# Output shot record as judiVector
Expand All @@ -45,7 +45,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Geometry, srcDa
rec_coords = setup_grid(recGeometry, modelPy.shape)

# Devito call
qOut = wrapcall_data(ac."adjoint_rec", modelPy, src_coords, rec_coords, dIn, space_order=options.space_order)
qOut = wrapcall_data(ac."adjoint_rec", modelPy, src_coords, rec_coords, dIn, space_order=options.space_order, f0=options.f0)
qOut = time_resample(qOut, dtComp, srcGeometry)

# Output adjoint data as judiVector
Expand All @@ -63,7 +63,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Geometry, srcDa
src_coords = setup_grid(srcGeometry, modelPy.shape)

# Devito call
u = wrapcall_function(ac."forward_no_rec", modelPy, src_coords, qIn, space_order=options.space_order)
u = wrapcall_function(ac."forward_no_rec", modelPy, src_coords, qIn, space_order=options.space_order, f0=options.f0)

# Output forward wavefield as judiWavefield
return judiWavefield(Info(prod(modelPy.shape), 1, size(u, 1)), dtComp, u)
Expand All @@ -80,7 +80,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Nothing, srcDat
rec_coords = setup_grid(recGeometry, modelPy.shape)

# Devito call
v = wrapcall_function(ac."adjoint_no_rec", modelPy, rec_coords, dIn, space_order=options.space_order)
v = wrapcall_function(ac."adjoint_no_rec", modelPy, rec_coords, dIn, space_order=options.space_order, f0=options.f0)

# Output adjoint wavefield as judiWavefield
return judiWavefield(Info(prod(modelPy.shape), 1, size(v, 1)), dtComp, v)
Expand All @@ -96,7 +96,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Nothing, srcDat
rec_coords = setup_grid(recGeometry, modelPy.shape)

# Devito call
dOut = wrapcall_data(ac."forward_wf_src", modelPy, srcData, rec_coords, space_order=options.space_order)
dOut = wrapcall_data(ac."forward_wf_src", modelPy, srcData, rec_coords, space_order=options.space_order, f0=options.f0)
dOut = time_resample(dOut, dtComp, recGeometry)

return judiVector{Float32, Array{Float32, 2}}("F*u", prod(size(dOut)), 1, 1, recGeometry, [dOut])
Expand All @@ -112,8 +112,8 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Geometry, srcDa
src_coords = setup_grid(srcGeometry, modelPy.shape)

# Devito call
qOut = wrapcall_data(ac."adjoint_wf_src", modelPy, recData, src_coords, space_order=options.space_order)
qOut = time_resample(qOut ,dtComp, srcGeometry)
qOut = wrapcall_data(ac."adjoint_wf_src", modelPy, recData, src_coords, space_order=options.space_order, f0=options.f0)
qOut = time_resample(qOut, dtComp, srcGeometry)

# Output adjoint data as judiVector
return judiVector{Float32, Array{Float32, 2}}("F'*d", prod(size(qOut)), 1, 1, srcGeometry, [qOut])
Expand All @@ -126,7 +126,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Nothing, srcDat
dtComp = convert(Float32, modelPy."critical_dt")

# Devito call
u = wrapcall_function(ac."forward_wf_src_norec", modelPy, srcData, space_order=options.space_order)
u = wrapcall_function(ac."forward_wf_src_norec", modelPy, srcData, space_order=options.space_order, f0=options.f0)

# Output forward wavefield as judiWavefield
return judiWavefield(Info(prod(modelPy.shape), 1, size(u, 1)), dtComp, u)
Expand All @@ -139,7 +139,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Nothing, srcDat
dtComp = convert(Float32, modelPy."critical_dt")

# Devito call
v = wrapcall_function(ac."adjoint_wf_src_norec", modelPy, recData, space_order=options.space_order)
v = wrapcall_function(ac."adjoint_wf_src_norec", modelPy, recData, space_order=options.space_order, f0=options.f0)

# Output adjoint wavefield as judiWavefield
return judiWavefield(Info(prod(modelPy.shape), 1, size(v, 1)), dtComp, v)
Expand All @@ -160,7 +160,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Geometry, srcDa

# Devito call
dOut = wrapcall_data(ac."born_rec", modelPy, src_coords, qIn, rec_coords,
space_order=options.space_order, isic=options.isic)
space_order=options.space_order, isic=options.isic, f0=options.f0)
dOut = time_resample(dOut, dtComp, recGeometry)

# Output linearized shot records as judiVector
Expand Down Expand Up @@ -188,8 +188,8 @@ function devito_interface(modelPy::PyCall.PyObject, srcGeometry::Geometry, srcDa
grad = wrapcall_function(ac."J_adjoint", modelPy,
src_coords, qIn, rec_coords, dIn, t_sub=options.subsampling_factor,
space_order=options.space_order, checkpointing=options.optimal_checkpointing,
freq_list=freqs, isic=options.isic,
dft_sub=options.dft_subsampling_factor[1])
freq_list=freqs, isic=options.isic, is_residual=true,
dft_sub=options.dft_subsampling_factor[1], f0=options.f0)

# Remove PML and return gradient as Array
grad = remove_padding(grad, modelPy.padsizes; true_adjoint=options.sum_padding)
Expand All @@ -211,8 +211,9 @@ function devito_interface(modelPy::PyCall.PyObject, srcData::Array, recGeometry:

# Devito call
dOut = wrapcall_data(ac."forward_rec_w", modelPy, weights,
qIn, rec_coords, space_order=options.space_order)
qIn, rec_coords, space_order=options.space_order, f0=options.f0)
dOut = time_resample(dOut, dtComp, recGeometry)

# Output shot record as judiVector
return judiVector{Float32, Array{Float32, 2}}("F*w", prod(size(dOut)), 1, 1, recGeometry, [dOut])
end
Expand All @@ -230,7 +231,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcData::Array, recGeometry:

# Devito call
wOut = wrapcall_function(ac."adjoint_w", modelPy, rec_coords, dIn,
qIn, space_order=options.space_order)
qIn, space_order=options.space_order, f0=options.f0)

# Output adjoint data as judiVector
wOut = remove_padding(wOut, modelPy.padsizes; true_adjoint=false)
Expand All @@ -253,7 +254,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcData::Array, recGeometry:

# Devito call
dOut = wrapcall_data(ac."born_rec_w", modelPy, weights, qIn, rec_coords,
space_order=options.space_order, isic=options.isic)
space_order=options.space_order, isic=options.isic, f0=options.f0)
dOut = time_resample(dOut, dtComp, recGeometry)

# Output linearized shot records as judiVector
Expand All @@ -266,7 +267,7 @@ function devito_interface(modelPy::PyCall.PyObject, srcData::Array, recGeometry:
end

# Adjoint Jacobian of extended source modeling: dm = J'*d_lin
function devito_interface(modelPy::PyCall.PyObject, srcData::Array, recGeometry::Geometry, recData::Array, weights:: Array, dm::Nothing, options::Options)
function devito_interface(modelPy::PyCall.PyObject, srcData::Array, recGeometry::Geometry, recData::Array, weights::Array, dm::Nothing, options::Options)

# Interpolate input data to computational grid
dtComp = convert(Float32, modelPy."critical_dt")
Expand All @@ -279,8 +280,8 @@ function devito_interface(modelPy::PyCall.PyObject, srcData::Array, recGeometry:
grad = wrapcall_function(ac."J_adjoint", modelPy,
nothing, qIn, rec_coords, dIn, t_sub=options.subsampling_factor,
space_order=options.space_order, checkpointing=options.optimal_checkpointing,
freq_list=freqs, isic=options.isic, ws=weights,
dft_sub=options.dft_subsampling_factor[1])
freq_list=freqs, isic=options.isic, ws=weights, is_residual=true,
dft_sub=options.dft_subsampling_factor[1], f0=options.f0)
# Remove PML and return gradient as Array
grad = remove_padding(grad, modelPy.padsizes; true_adjoint=options.sum_padding)
return PhysicalParameter(grad, modelPy.spacing, modelPy.origin)
Expand Down
2 changes: 1 addition & 1 deletion src/TimeModeling/Modeling/twri_objective_serial.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ function twri_objective(model_full::Model, source::judiVector, dObs::judiVector,
t_sub=options.subsampling_factor, space_order=options.space_order,
grad=optionswri.params, grad_corr=optionswri.grad_corr, eps=optionswri.eps,
alpha_op=optionswri.comp_alpha, w_fun=optionswri.weight_fun,
freq_list=freqs, wfilt=wfilt)
freq_list=freqs, wfilt=wfilt, f0=options.f0)

if (optionswri.params in [:m, :all])
gradm = remove_padding(gradm, modelPy.padsizes; true_adjoint=options.sum_padding)
Expand Down
10 changes: 7 additions & 3 deletions src/TimeModeling/Types/ModelStructure.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -299,10 +299,10 @@ end

function Model(n::IntTuple, d::RealTuple, o::RealTuple, m;
epsilon=nothing, delta=nothing, theta=nothing,
phi=nothing, rho=nothing, nb=40)
phi=nothing, rho=nothing, qp=nothing, nb=40)

params = Dict(:m => PhysicalParameter(Float32.(m), d, o))
for (name, val)=zip([:rho, :epsilon, :delta, :theta, :phi], [rho, epsilon, delta, theta, phi])
for (name, val)=zip([:rho, :qp, :epsilon, :delta, :theta, :phi], [rho, qp, epsilon, delta, theta, phi])
~isnothing(val) && (params[name] = PhysicalParameter(Float32.(val), n, d, o))
end

Expand All @@ -313,6 +313,10 @@ function Model(n::IntTuple, d::RealTuple, o::RealTuple, m, rho; nb=40)
return Model(n, d, o, m; rho=rho, nb=nb)
end

function Model(n::IntTuple, d::RealTuple, o::RealTuple, m, rho, qp; nb=40)
return Model(n, d, o, m; rho=rho, qp=qp, nb=nb)
end

get_dt(m::Model; dt=nothing) = calculate_dt(m; dt=dt)
getindex(m::Model, sym::Symbol) = m.params[sym]

Expand All @@ -330,4 +334,4 @@ function Base.getproperty(obj::Model, sym::Symbol)
end

similar(x::PhysicalParameter{vDT}, m::Model) where {vDT} = PhysicalParameter(m.n, m.d, m.o; vDT=vDT)
similar(x::Array, m::Model) where {vDT} = similar(x, m.n)
similar(x::Array, m::Model) where {vDT} = similar(x, m.n)
Loading

1 comment on commit e2307a4

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/58829

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v2.6.9 -m "<description of version>" e2307a443f75ffde1451a1916410559e512eb249
git push origin v2.6.9

Please sign in to comment.