diff --git a/.github/workflows/ci-judi.yml b/.github/workflows/ci-judi.yml index 073802aca..9158927c0 100644 --- a/.github/workflows/ci-judi.yml +++ b/.github/workflows/ci-judi.yml @@ -19,7 +19,7 @@ jobs: name: JUDI base on Julia ${{ matrix.version }} - ${{ matrix.os }} runs-on: ${{ matrix.os }} env: - DEVITO_ARCH: gcc-11 + DEVITO_ARCH: gcc-12 DEVITO_LANGUAGE: "openmp" OMP_NUM_THREADS: 4 GROUP: "JUDI" diff --git a/.github/workflows/ci-op.yml b/.github/workflows/ci-op.yml index 847190ff7..a393e4995 100644 --- a/.github/workflows/ci-op.yml +++ b/.github/workflows/ci-op.yml @@ -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: ${{ matrix.cc }} DEVITO_LANGUAGE: "openmp" DEVITO_LOGGING: "INFO" OMP_NUM_THREADS: ${{ matrix.omp }} @@ -33,32 +33,38 @@ jobs: op: ["ISO_OP", "ISO_OP_FS", "TTI_OP", "TTI_OP_FS", "BASICS"] version: ['1'] omp: [2] + cc: ['gcc-11'] include: - os: macos-13 version: '1.6' op: "ISO_OP" omp: 1 + cc: gcc-13 - os: macos-13 version: '1.8' op: "ISO_OP" omp: 1 + cc: gcc-13 - os: macos-13 version: '1.9' op: "ISO_OP" omp: 1 + cc: gcc-13 - os: ubuntu-latest version: '1.9' op: "VISCO_AC_OP" omp: 2 + cc: gcc-11 - os: ubuntu-latest version: '1.10' op: "ISO_OP" omp: 2 + cc: gcc-11 steps: - name: Checkout JUDI diff --git a/Project.toml b/Project.toml index 06e036797..d239b0046 100644 --- a/Project.toml +++ b/Project.toml @@ -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" diff --git a/examples/scripts/modeling_basic_elastic.jl b/examples/scripts/modeling_basic_elastic.jl index 454e67561..c252f9606 100644 --- a/examples/scripts/modeling_basic_elastic.jl +++ b/examples/scripts/modeling_basic_elastic.jl @@ -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] @@ -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] @@ -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 diff --git a/src/pysource/geom_utils.py b/src/pysource/geom_utils.py index 5f8e49c06..0dccc325b 100644 --- a/src/pysource/geom_utils.py +++ b/src/pysource/geom_utils.py @@ -7,7 +7,7 @@ try: from recipes.utils import mirror_source except ImportError: - def mirror_source(src): + def mirror_source(model, src): return src @@ -20,12 +20,12 @@ def src_rec(model, u, src_coords=None, rec_coords=None, wavelet=None, nt=None): src = wavelet else: src = PointSource(name="src%s" % namef, grid=model.grid, ntime=nt, - coordinates=src_coords, interpolation='sinc', r=3) + coordinates=src_coords) src.data[:] = wavelet.view(np.ndarray) if wavelet is not None else 0. rcv = None if rec_coords is not None: rcv = Receiver(name="rcv%s" % namef, grid=model.grid, ntime=nt, - coordinates=rec_coords, interpolation='sinc', r=3) + coordinates=rec_coords) return src, rcv @@ -62,14 +62,20 @@ 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) + + geom_expr = [] + # Source 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 @@ -77,13 +83,14 @@ def geom_expr(model, u, src_coords=None, rec_coords=None, wavelet=None, fw=True, if model.fs: # Free surface src_eq = mirror_source(model, src_eq) - geom_expr += src_eq + + geom_expr += [src_eq] # Setup adjoint wavefield sampling at source locations if rcv is not None: if model.is_elastic: - rec_expr = u[1].trace() + geom_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 + geom_expr += rcv.interpolate(expr=rec_expr) + return geom_expr diff --git a/src/pysource/models.py b/src/pysource/models.py index 9bf958fa9..81f3f069b 100644 --- a/src/pysource/models.py +++ b/src/pysource/models.py @@ -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'] @@ -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 @@ -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 @@ -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 @@ -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 @@ -313,12 +290,12 @@ def zero_thomsen(self): def _gen_phys_param(self, field, name, space_order, is_param=False, default_value=0, avg_mode='arithmetic'): """ - Create symbolic object an initiliaze its data + Create symbolic object an initialize its data """ 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: @@ -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)) @@ -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: @@ -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] @@ -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): """ @@ -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 @@ -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 @@ -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) diff --git a/src/pysource/operators.py b/src/pysource/operators.py index 16b464b43..cd99036ca 100644 --- a/src/pysource/operators.py +++ b/src/pysource/operators.py @@ -99,8 +99,8 @@ def forward_op(p_params, tti, visco, elas, space_order, fw, spacing, save, t_sub u = wavefield(model, space_order, save=save, nt=nt, t_sub=t_sub, fw=fw) # Setup source and receiver - g_expr = geom_expr(model, u, src_coords=scords, nt=nt, - rec_coords=rcords, wavelet=wavelet, fw=fw) + gexpr = geom_expr(model, u, src_coords=scords, nt=nt, + rec_coords=rcords, wavelet=wavelet, fw=fw) # Expression for saving wavefield if time subsampling is used eq_save = save_subsampled(model, u, nt, t_sub, space_order=space_order) @@ -126,7 +126,7 @@ def forward_op(p_params, tti, visco, elas, space_order, fw, spacing, save, t_sub # Create operator and run subs = model.spacing_map pname = "forward" if fw else "adjoint" - op = Operator(pde + wrec + nv_t + dft + g_expr + extra + eq_save + nv_s + Ieq, + op = Operator(pde + wrec + nv_t + dft + gexpr + extra + eq_save + nv_s + Ieq, subs=subs, name=pname+name(model), opt=opt_op(model)) op.cfunction @@ -158,9 +158,9 @@ def born_op(p_params, tti, visco, elas, space_order, fw, spacing, save, pt_src, ul = wavefield(model, space_order, name="l", fw=fw) # Setup source and receiver - g_expr = geom_expr(model, u, rec_coords=rcords if nlind else None, - src_coords=scords, wavelet=wavelet, fw=fw) - g_exprl = geom_expr(model, ul, rec_coords=rcords, nt=nt, fw=fw) + gexpr = geom_expr(model, u, rec_coords=rcords if nlind else None, + src_coords=scords, wavelet=wavelet, fw=fw) + gexprl = geom_expr(model, ul, rec_coords=rcords, nt=nt, fw=fw) # Expression for saving wavefield if time subsampling is used eq_save = save_subsampled(model, u, nt, t_sub, space_order=space_order) @@ -183,7 +183,7 @@ def born_op(p_params, tti, visco, elas, space_order, fw, spacing, save, pt_src, # Create operator and run subs = model.spacing_map - op = Operator(pde + g_expr + extra + g_exprl + pdel + extral + dft + eq_save + Ieq, + op = Operator(pde + gexpr + extra + gexprl + pdel + extral + dft + eq_save + Ieq, subs=subs, name="born"+name(model), opt=opt_op(model)) op.cfunction @@ -211,7 +211,7 @@ def adjoint_born_op(p_params, tti, visco, elas, space_order, fw, spacing, pt_rec dft=nfreq > 0, t_sub=t_sub, fw=fw) # Setup source and receiver - r_expr = geom_expr(model, v, src_coords=rcords, wavelet=residual, fw=not fw) + gexpr = geom_expr(model, v, src_coords=rcords, wavelet=residual, fw=not fw) # Set up PDE expression and rearrange pde, extra = wave_kernel(model, v, fw=False, f0=Constant('f0')) @@ -226,7 +226,7 @@ def adjoint_born_op(p_params, tti, visco, elas, space_order, fw, spacing, pt_rec # Create operator and run subs = model.spacing_map - op = Operator(pde + r_expr + extra + g_expr + Ieq, subs=subs, + op = Operator(pde + gexpr + extra + g_expr + Ieq, subs=subs, name="gradient"+name(model), opt=opt_op(model)) try: diff --git a/src/pysource/propagators.py b/src/pysource/propagators.py index 1cdea79ff..c6ebc1978 100644 --- a/src/pysource/propagators.py +++ b/src/pysource/propagators.py @@ -253,7 +253,7 @@ def forward_grad(model, src_coords, rcv_coords, wavelet, v, pde, extra = wave_kernel(model, u, q=q, f0=f0) # Setup source and receiver - rexpr = geom_expr(model, u, src_coords=src_coords, nt=nt, + gexpr = geom_expr(model, u, src_coords=src_coords, nt=nt, rec_coords=rcv_coords, wavelet=wavelet) _, rcv = src_rec(model, u, src_coords, rcv_coords, wavelet, nt) @@ -263,7 +263,7 @@ def forward_grad(model, src_coords, rcv_coords, wavelet, v, # Create operator and run subs = model.spacing_map - op = Operator(pde + rexpr + extra + g_expr, + op = Operator(pde + gexpr + extra + g_expr, subs=subs, name="forward_grad", opt=opt_op(model))