Skip to content

Commit

Permalink
need vp
Browse files Browse the repository at this point in the history
  • Loading branch information
@mloubout
mloubout committed Sep 30, 2024
1 parent af00642 commit da70e62
Show file tree
Hide file tree
Showing 7 changed files with 80 additions and 158 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/ci-op.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ jobs:
name: ${{ matrix.op }} on Julia ${{ matrix.version }} - ${{ matrix.os }}
runs-on: ${{ matrix.os }}
env:
DEVITO_ARCH: gcc-11
DEVITO_ARCH: gcc-12
DEVITO_LANGUAGE: "openmp"
DEVITO_LOGGING: "INFO"
OMP_NUM_THREADS: ${{ matrix.omp }}
Expand Down
1 change: 1 addition & 0 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
JOLI = "bb331ad6-a1cf-11e9-23da-9bcb53c69f6f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Requires = "ae029012-a4dd-5104-9daa-d747884805df"
Expand Down
6 changes: 3 additions & 3 deletions examples/scripts/modeling_basic_elastic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ nxrec = 120
nyrec = 100
xrec = range(50f0, stop=1150f0, length=nxrec)
yrec = 0f0
zrec = range(0f0, stop=0f0, length=nxrec)
zrec = range(10f0, stop=10f0, length=nxrec)

# receiver sampling and recording time
timeR = 1500f0 # receiver recording time [ms]
Expand All @@ -41,7 +41,7 @@ recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)
# Set up source geometry (cell array with source locations for each shot)
xsrc = 600f0
ysrc = 0f0
zsrc = 0f0
zsrc = 10f0

# source sampling and number of time steps
timeS = 1500f0 # source length in [ms]
Expand All @@ -57,7 +57,7 @@ wavelet = ricker_wavelet(timeS, dtS, f0)
###################################################################################################

# Setup operators
F = judiModeling(model, srcGeometry, recGeometry)
F = judiModeling(model, srcGeometry, recGeometry; options=Options(space_order=8, free_surface=true))
q = judiVector(srcGeometry, wavelet)

# Nonlinear modeling
Expand Down
23 changes: 15 additions & 8 deletions src/pysource/geom_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
try:
from recipes.utils import mirror_source
except ImportError:
def mirror_source(src):
def mirror_source(model, src):
return src


Expand Down Expand Up @@ -62,28 +62,35 @@ def geom_expr(model, u, src_coords=None, rec_coords=None, wavelet=None, fw=True,
if not model.is_elastic:
m = model.m * irho
dt = model.grid.time_dim.spacing
geom_expr = []
src, rcv = src_rec(model, u, src_coords, rec_coords, wavelet, nt)
model.__init_abox__(src, rcv, fw=fw)

# Source
src_eq = []
if src is not None:
# Elastic inject into diagonal of stress
if model.is_elastic:
for ud in as_tuple(u)[1].diagonal():
geom_expr += src.inject(field=ud.forward, expr=src*dt/irho)
c = 1 / model.grid.dim
src_eq = src.inject(field=as_tuple(u)[1].forward.diagonal(),
expr=c*src*dt/irho)
if model.fs:
# Free surface
src_eq = mirror_source(model, src_eq)
else:
# Acoustic inject into pressure
u_n = as_tuple(u)[0].forward if fw else as_tuple(u)[0].backward
src_eq = src.inject(field=u_n, expr=src*dt**2/m)
if model.fs:
# Free surface
src_eq = mirror_source(model, src_eq)
geom_expr += src_eq

# Setup adjoint wavefield sampling at source locations
adj_rcv = []
if rcv is not None:
if model.is_elastic:
rec_expr = u[1].trace()
rec_expr = u[1].trace() / model.grid.dim
else:
rec_expr = u[0] if model.is_tti else u
adj_rcv = rcv.interpolate(expr=rec_expr)
geom_expr += adj_rcv
return geom_expr

return src_eq, adj_rcv
182 changes: 48 additions & 134 deletions src/pysource/models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,56 +9,18 @@

from devito.tools import as_tuple, memoized_func

try:
from devitopro import * # noqa
from devitopro.subdomains.abox import ABox, ABoxFunction
from devitopro.data import Float16
AboxBase = ABox
except ImportError:
ABox = None
AboxBase = object
ABoxFunction = object
Float16 = lambda *ar, **kw: np.float32
# try:
# from devitopro import * # noqa
# from devitopro.subdomains.abox import ABox, ABoxFunction
# from devitopro.data import Float16
# AboxBase = ABox
# except ImportError:
ABox = None
AboxBase = object
ABoxFunction = object
Float16 = lambda *ar, **kw: np.float32


class SlownessABoxFunction(ABoxFunction):

_physical_params = ('m',)

def vp_max(self, rdim, **kwargs):
vp = kwargs.get('vp', self.vp)
eps = kwargs.get('eps', self.eps)

vpi = vp.data.min(axis=rdim)**(-.5)
# Thomsen correction
if eps is not None:
epsi = eps.data.max(axis=rdim)
vpi._local[:] *= np.sqrt(1. + 2.*epsi._local[:])

return vpi


class LameABoxFunction(ABoxFunction):

_physical_params = ('mu', 'lam', 'b')

def vp_max(self, rdim, **kwargs):
lam = kwargs.get('lam', self.lam)
mu = kwargs.get('mu', self.mu)
b = kwargs.get('b', self.b)

return np.sqrt(((lam.data + 2 * mu.data) * b.data).max(axis=rdim))


class JUDIAbox(AboxBase):

def __init__(self, *args, subdomains=None, name=None, **params):
if 'mu' in params:
self._afunc = LameABoxFunction
else:
self._afunc = SlownessABoxFunction

super().__init__(*args, subdomains=subdomains, name=name, **params)
_dtypes = {'params': 'f32'}


__all__ = ['Model']
Expand Down Expand Up @@ -196,7 +158,7 @@ def __init__(self, origin, spacing, shape, space_order=8, nbl=40, dtype=np.float
self.shape = tuple(shape)
self.nbl = int(nbl)
self.origin = tuple([dtype(o) for o in origin])
abc_type = "mask" if (qp is not None or mu is not None) else "damp"
abc_type = "mask" if qp is not None else "damp"
self.fs = fs
self._abox = None
# Topology in case. Always decompose only in x or y
Expand Down Expand Up @@ -226,13 +188,17 @@ def __init__(self, origin, spacing, shape, space_order=8, nbl=40, dtype=np.float
# Seismic fields and properties
self.scale = 1
self._space_order = space_order

# Create square slowness of the wave as symbol `m`
if m is not None:
self._m = self._gen_phys_param(m, 'm', space_order)
vp_vals = m**(-.5)
self.m = self._gen_phys_param(m, 'm', space_order)

# density
self._init_density(rho, b, space_order)

# Perturbation for linearized modeling
self._dm = self._gen_phys_param(dm, 'dm', space_order)
self.dm = self._gen_phys_param(dm, 'dm', space_order)

# Model type
self._is_viscoacoustic = qp is not None
Expand All @@ -257,11 +223,19 @@ def __init__(self, origin, spacing, shape, space_order=8, nbl=40, dtype=np.float
# Additional parameter fields for elastic
if self._is_elastic:
self.lam = self._gen_phys_param(lam, 'lam', space_order, is_param=True)
mu[np.where(mu == 0)] = 1e-12
self.mu = self._gen_phys_param(mu, 'mu', space_order, is_param=True,
avg_mode='harmonic')
b = b if b is not None else 1 / rho
vp_vals = ((lam + 2 * mu) * b)**(.5)

# User provided dt
self._dt = kwargs.get('dt')

# Need vp for Abox
if ABox is not None:
self._vp = self._gen_phys_param(vp_vals, '_vp', space_order)

def _init_density(self, rho, b, so):
"""
Initialize density parameter. Depending on variance in density
Expand Down Expand Up @@ -297,6 +271,9 @@ def physical_params(self, **kwargs):
if not kwargs.get('born', False):
params.pop('dm', None)

# Remove "build" _x params
params.pop('_vp', None)

return params

@property
Expand All @@ -318,7 +295,7 @@ def _gen_phys_param(self, field, name, space_order, is_param=False,
if field is None:
return default_value
if isinstance(field, np.ndarray):
if field.dtype == np.float16:
if _dtypes['params'] == 'f16' or field.dtype == np.float16:
_min = np.amin(field)
_max = np.amax(field)
if _max == _min:
Expand All @@ -329,7 +306,7 @@ def _gen_phys_param(self, field, name, space_order, is_param=False,

function = Function(name=name, grid=self.grid, space_order=space_order,
parameter=is_param, avg_mode=avg_mode, dtype=dtype)
pad = self.padsizes if field.shape == self.shape else 0
pad = 0 if field.shape == function.shape else self.padsizes
initialize_function(function, field, pad)
else:
function = Constant(name=name, value=np.amin(field))
Expand All @@ -346,9 +323,12 @@ def physical_parameters(self):
param = getattr(self, p)
# Get dtype
comp = getattr(param, '_compression', param.dtype)
if isinstance(comp, Float16):
dtype = comp
else:
try:
if isinstance(comp, Float16):
dtype = comp
else:
dtype = param.dtype
except TypeError:
dtype = param.dtype
# Add to list
if param.is_Constant:
Expand Down Expand Up @@ -457,7 +437,7 @@ def _cfl_coeff(self):
so = max(self.space_order // 2, 2)
coeffs = fd_w(1, range(-so, so), .5)
c_fd = sum(np.abs(coeffs[-1][-1])) / 2
return .9 * np.sqrt(self.dim) / self.dim / c_fd
return .95 * np.sqrt(self.dim) / self.dim / c_fd
a1 = 4 # 2nd order in time
so = max(self.space_order // 2, 4)
coeffs = fd_w(2, range(-so, so), 0)[-1][-1]
Expand Down Expand Up @@ -489,71 +469,6 @@ def critical_dt(self):
return self.dtype("%.3e" % self.dt)
return dt

@property
def dm(self):
"""
Model perturbation for linearized modeling
"""
return self._dm

@dm.setter
def dm(self, dm):
"""
Set a new model perturbation.
Parameters
----------
dm : float or array
New model perturbation
"""
# Update the square slowness according to new value
if isinstance(dm, np.ndarray):
if not isinstance(self._dm, Function):
self._dm = self._gen_phys_param(dm, 'dm', self.space_order)
elif dm.shape == self.shape:
initialize_function(self._dm, dm, self.padsizes)
elif dm.shape == self.dm.shape:
self.dm.data[:] = dm[:]
else:
raise ValueError("Incorrect input size %s for model of size" % dm.shape +
" %s without or %s with padding" % (self.shape,
self.dm.shape))
else:
try:
self._dm.data = dm
except AttributeError:
self._dm = dm

@property
def m(self):
"""
Function holding the squared slowness in s^2/km^2.
"""
return self._m

@m.setter
def m(self, m):
"""
Set a new squared slowness model.
Parameters
----------
m : float or array
New squared slowness in s^2/km^2.
"""
# Update the square slowness according to new value
if isinstance(m, np.ndarray):
if m.shape == self.m.shape:
self.m.data[:] = m[:]
elif m.shape == self.shape:
initialize_function(self._m, m, self.padsizes)
else:
raise ValueError("Incorrect input size %s for model of size" % m.shape +
" %s without or %s with padding" % (self.shape,
self.m.shape))
else:
self._m.data = m

@property
def vp(self):
"""
Expand All @@ -575,8 +490,9 @@ def abox(self, src, rec, fw=True):
return {}
if not fw:
src, rec = rec, src
abox = JUDIAbox(self.space_order, src=src, rcv=rec, **self.physical_params())
return {'abox': abox._abox_func}
eps = getattr(self, 'epsilon', None)
abox = ABox(src, rec, self._vp, self.space_order, eps=eps)
return {'abox': abox}

def __init_abox__(self, src, rec, fw=True):
return
Expand Down Expand Up @@ -622,15 +538,16 @@ def __init__(self, tti, visco, elastic, spacing, fs, space_order, p_params):
self.damp = Function(name='damp', grid=self.grid, space_order=0)
_physical_parameters = ['damp']
for p, dt in p_params.items():
if _dtypes['params'] == 'f16':
dt = np.float16
if p.endswith('_const'):
name = p.split('_')[0]
setattr(self, name, Constant(name=name, value=1, dtype=dt))
else:
pn = '_%s' % p if p in ['m', 'dm'] else p
avgmode = 'harmonic' if p == 'mu' else 'arithmetic'
setattr(self, pn, Function(name=p, grid=self.grid, is_param=True,
space_order=space_order, dtype=dt,
avg_mode=avgmode))
setattr(self, p, Function(name=p, grid=self.grid, is_param=True,
space_order=space_order, dtype=dt,
avg_mode=avgmode))
_physical_parameters.append(p)
if 'irho' not in p_params and 'irho_const' not in p_params:
self.irho = 1 if 'rho' not in p_params else 1 / self.rho
Expand Down Expand Up @@ -663,8 +580,5 @@ def __init_abox__(self, src, rec, fw=True):

if not fw:
src, rec = rec, src
try:
self._abox = JUDIAbox(self.space_order, src=src, rcv=rec,
**self.physical_params())
except AttributeError:
return
eps = getattr(self, 'epsilon', None)
self._abox = ABox(src, rec, self._vp, self.space_order, eps=eps)
Loading

0 comments on commit da70e62

Please sign in to comment.