Refine CFL condition in seismic

devito/types/basic.py

@@ -16,7 +16,7 @@
 from devito.parameters import configuration
 from devito.symbolics import aligned_indices
 from devito.tools import (Pickable, ctypes_to_cstr, dtype_to_cstr, dtype_to_ctype,
-                          frozendict)
+                          frozendict, memoized_meth)
 from devito.types.args import ArgProvider
 from devito.types.caching import Cached
 from devito.types.utils import DimensionTuple
@@ -729,15 +729,23 @@ def dimensions(self):
     def _eval_deriv(self):
         return self
 
-    @cached_property
+    @property
     def _is_on_grid(self):
         """
-        Check whether the object is on the grid or need averaging.
+        Check whether the object is on the grid and requires averaging.
         For example, if the original non-staggered function is f(x)
         then f(x) is on the grid and f(x + h_x/2) is off the grid.
         """
+        return self._check_indices(inds=self.indices)
+
+    @memoized_meth
+    def _check_indices(self, inds=None):
+        """
+        Check if the function indices are aligned with the dimensions.
+        """
+        inds = inds or self.indices
         return all([aligned_indices(i, j, d.spacing) for i, j, d in
-                    zip(self.indices, self.indices_ref, self.dimensions)])
+                    zip(inds, self.indices_ref, self.dimensions)])
 
     @property
     def evaluate(self):

examples/seismic/acoustic/acoustic_example.py

@@ -1,17 +1,18 @@
 import numpy as np
+import pytest
 
 from devito.logger import info
-from devito import Constant, Function, smooth
+from devito import Constant, Function, smooth, norm
 from examples.seismic.acoustic import AcousticWaveSolver
 from examples.seismic import demo_model, setup_geometry, seismic_args
 
 
 def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
                    tn=500., kernel='OT2', space_order=4, nbl=10,
-                   preset='layers-isotropic', **kwargs):
+                   preset='layers-isotropic', fs=False, **kwargs):
     model = demo_model(preset, space_order=space_order, shape=shape, nbl=nbl,
                        dtype=kwargs.pop('dtype', np.float32), spacing=spacing,
-                       **kwargs)
+                       fs=fs, **kwargs)
 
     # Source and receiver geometries
     geometry = setup_geometry(model, tn)
@@ -23,11 +24,11 @@ def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
 
 
 def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
-        space_order=4, kernel='OT2', nbl=40, full_run=False,
+        space_order=4, kernel='OT2', nbl=40, full_run=False, fs=False,
         autotune=False, preset='layers-isotropic', checkpointing=False, **kwargs):
 
     solver = acoustic_setup(shape=shape, spacing=spacing, nbl=nbl, tn=tn,
-                            space_order=space_order, kernel=kernel,
+                            space_order=space_order, kernel=kernel, fs=fs,
                             preset=preset, **kwargs)
 
     info("Applying Forward")
@@ -62,6 +63,20 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
     return summary.gflopss, summary.oi, summary.timings, [rec, u.data]
 
 
+@pytest.mark.parametrize('ndim', [1, 2, 3])
+def test_isoacoustic_stability(ndim):
+    shape = tuple([11]*ndim)
+    spacing = tuple([20]*ndim)
+    _, _, _, [rec, _] = run(shape=shape, spacing=spacing, tn=20000.0, nbl=0)
+    assert np.isfinite(norm(rec))
+
+
+@pytest.mark.parametrize('fs, normrec', [(True, 369.955), (False, 459.1678)])
+def test_isoacoustic(fs, normrec):
+    _, _, _, [rec, _] = run(fs=fs)
+    assert np.isclose(norm(rec), normrec, rtol=1e-3, atol=0)
+
+
 if __name__ == "__main__":
     description = ("Example script for a set of acoustic operators.")
     parser = seismic_args(description)
@@ -76,7 +91,7 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
     tn = args.tn if args.tn > 0 else (750. if ndim < 3 else 1250.)
 
     preset = 'constant-isotropic' if args.constant else 'layers-isotropic'
-    run(shape=shape, spacing=spacing, nbl=args.nbl, tn=tn,
+    run(shape=shape, spacing=spacing, nbl=args.nbl, tn=tn, fs=args.fs,
         space_order=args.space_order, preset=preset, kernel=args.kernel,
         autotune=args.autotune, opt=args.opt, full_run=args.full,
         checkpointing=args.checkpointing)

examples/seismic/acoustic/operators.py

@@ -17,7 +17,7 @@ def freesurface(model, eq):
     eq : Eq
         Time-stepping stencil (time update) to mirror at the freesurface.
     """
-    lhs, rhs = eq.evaluate
+    lhs, rhs = eq.evaluate.args
     # Get vertical dimension and corresponding subdimension
     zfs = model.grid.subdomains['fsdomain'].dimensions[-1]
     z = zfs.parent

examples/seismic/acoustic/wavesolver.py

@@ -36,9 +36,9 @@ def __init__(self, model, geometry, kernel='OT2', space_order=4, **kwargs):
         self.kernel = kernel
 
         # Time step can be \sqrt{3}=1.73 bigger with 4th order
-        self.dt = self.model.critical_dt
         if self.kernel == 'OT4':
-            self.dt = model.dtype(self.dt*1.73)
+            self.model.dt_scale = 1.73
+        self.dt = self.model.critical_dt
 
         # Cache compiler options
         self._kwargs = kwargs

examples/seismic/elastic/elastic_example.py

@@ -1,6 +1,8 @@
 import numpy as np
+import pytest
 
 from devito.logger import info
+from devito import norm
 from examples.seismic.elastic import ElasticWaveSolver
 from examples.seismic import demo_model, setup_geometry, seismic_args
 
@@ -35,9 +37,16 @@ def run(shape=(50, 50), spacing=(20.0, 20.0), tn=1000.0,
 
 def test_elastic():
     _, _, _, [rec1, rec2, v, tau] = run()
-    norm = lambda x: np.linalg.norm(x.data.reshape(-1))
-    assert np.isclose(norm(rec1), 20.59193, atol=1e-3, rtol=0)
-    assert np.isclose(norm(rec2), 0.671578, atol=1e-3, rtol=0)
+    assert np.isclose(norm(rec1), 19.25636, atol=1e-3, rtol=0)
+    assert np.isclose(norm(rec2), 0.627606, atol=1e-3, rtol=0)
+
+
+@pytest.mark.parametrize('ndim', [1, 2, 3])
+def test_elastic_stability(ndim):
+    shape = tuple([11]*ndim)
+    spacing = tuple([20]*ndim)
+    _, _, _, [rec1, rec2, v, tau] = run(shape=shape, spacing=spacing, tn=20000.0, nbl=0)
+    assert np.isfinite(norm(rec1))
 
 
 if __name__ == "__main__":

examples/seismic/model.py

@@ -1,6 +1,5 @@
 import numpy as np
-from sympy import sin, Abs, finite_diff_weights
-
+from sympy import sin, Abs, finite_diff_weights as fd_w
 
 from devito import (Grid, SubDomain, Function, Constant,
                     SubDimension, Eq, Inc, Operator, div)
@@ -77,7 +76,7 @@ def __init__(self, so):
 
     def define(self, dimensions):
         """
-        Definition of the top part of the domain for wrapped indices FS
+        Definition of the upper section of the domain for wrapped indices FS.
         """
 
         return {d: (d if not d == dimensions[-1] else ('left', self.size))
@@ -259,6 +258,11 @@ def __init__(self, origin, spacing, shape, space_order, vp, nbl=20, fs=False,
 
         # User provided dt
         self._dt = kwargs.get('dt')
+        # Some wave equation need a rescaled dt that can't be infered from the model
+        # parameters, such as isoacoustic OT4 that can use a dt sqrt(3) larger than
+        # isoacoustic OT2. This property should be set from a wavesolver or after model
+        # instanciation only via model.dt_scale = value.
+        self._dt_scale = 1
 
     def _initialize_physics(self, vp, space_order, **kwargs):
         """
@@ -299,31 +303,36 @@ def _max_vp(self):
             return np.sqrt(mmin(self.b) * (mmax(self.lam) + 2 * mmax(self.mu)))
 
     @property
-    def _scale(self):
+    def _thomsen_scale(self):
         # Update scale for tti
         if 'epsilon' in self._physical_parameters:
             return np.sqrt(1 + 2 * mmax(self.epsilon))
         return 1
 
+    @property
+    def dt_scale(self):
+        return self._dt_scale
+
+    @dt_scale.setter
+    def dt_scale(self, val):
+        self._dt_scale = val
+
     @property
     def _cfl_coeff(self):
         """
-        Courant number from the physics.
-        For the elastic case, we use the general `.85/sqrt(ndim)`.
-
-        For te acoustic case, we can get an accurate estimate from the spatial order.
-        One dan show that C = sqrt(a1/a2) where:
-        - a1 is the sum of absolute value of FD coefficients in time
-        - a2 is the sum of absolute value of FD coefficients in space
-        https://library.seg.org/doi/pdf/10.1190/1.1444605
+        Courant number from the physics and spatial discretization order.
+        The CFL coefficients are described in:
+        - https://doi.org/10.1137/0916052 for the elastic case
+        - https://library.seg.org/doi/pdf/10.1190/1.1444605 for the acoustic case
         """
         # Elasic coefficient (see e.g )
         if 'lam' in self._physical_parameters or 'vs' in self._physical_parameters:
-            return (.85 if self.grid.dim == 3 else .75) / np.sqrt(self.grid.dim)
+            coeffs = fd_w(1, range(-self.space_order//2+1, self.space_order//2+1), .5)
+            c_fd = sum(np.abs(coeffs[-1][-1])) / 2
+            return np.sqrt(self.dim) / self.dim / c_fd
         a1 = 4  # 2nd order in time
-        coeffs = finite_diff_weights(2, range(-self.space_order, self.space_order+1), 0)
-        a2 = float(self.grid.dim * sum(np.abs(coeffs[-1][-1])))
-        return np.sqrt(a1/a2)
+        coeffs = fd_w(2, range(-self.space_order, self.space_order+1), 0)[-1][-1]
+        return np.sqrt(a1/float(self.grid.dim * sum(np.abs(coeffs))))
 
     @property
     def critical_dt(self):
@@ -334,7 +343,8 @@ def critical_dt(self):
         #
         # The CFL condtion is then given by
         # dt <= coeff * h / (max(velocity))
-        dt = self._cfl_coeff * np.min(self.spacing) / (self._scale*self._max_vp)
+        dt = self._cfl_coeff * np.min(self.spacing) / (self._thomsen_scale*self._max_vp)
+        dt = self.dt_scale * dt
         if self._dt:
             assert self._dt < dt
             return self._dt

examples/seismic/skew_self_adjoint/wavesolver.py

@@ -46,7 +46,8 @@ def __init__(self, model, geometry, space_order=8, **kwargs):
         self.space_order = space_order
 
         # Time step is .6 time smaller due to Q
-        self.dt = .6*self.model.critical_dt
+        self.model.dt_scale = .6
+        self.dt = self.model.critical_dt
 
         # Cache compiler options
         self._kwargs = kwargs