model: make dt user side scalable with dt_scale

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/model.py

@@ -258,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):
         """
@@ -298,12 +303,20 @@ 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):
         """
@@ -330,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

examples/seismic/viscoelastic/wavesolver.py

@@ -29,7 +29,12 @@ def __init__(self, model, geometry, space_order=4, **kwargs):
         self.geometry = geometry
 
         self.space_order = space_order
-        self.dt = self.model.dtype(.9 * self.model.critical_dt)
+
+        # The viscoelastic equation requires a smaller dt than the standard
+        # elastic equation due to instability introduced by the viscosity.
+        self.model.dt_scale = .9
+        self.dt = self.model.critical_dt
+
         # Cache compiler options
         self._kwargs = kwargs