Merge pull request #231 from prabhuramachandran/dt_adapt

Simpler adaptive timestep

docs/source/design/equations.rst

@@ -447,10 +447,45 @@ equation in the list of equations to be passed to the application.  It is
 often easiest to look at the many existing equations in PySPH and learn the
 general patterns.
 
-If you wish to use adaptive time stepping, see the code
-:py:class:`pysph.sph.integrator.Integrator`. The integrator uses information
-from the arrays ``dt_cfl``, ``dt_force``, and ``dt_visc`` in each of the
-particle arrays to determine the most suitable time step.
+Adaptive timesteps
+--------------------
+
+There are a couple of ways to use adaptive timesteps. The first is to compute
+a required timestep directly per-particle in a particle array property called
+``dt_adapt``. The minimum value of this array across all particle arrays is
+used to set the timestep directly. This is the easiest way to set the adaptive
+timestep.
+
+If the ``dt_adapt`` parameter is not set one may also use standard velocity,
+force, and viscosity based parameters. The integrator uses information from
+the arrays ``dt_cfl``, ``dt_force``, and ``dt_visc`` in each of the particle
+arrays to determine the most suitable time step. This is done using the
+following approach. The minimum smoothing parameter ``h`` is found as
+``hmin``. Let the CFL number be given as ``cfl``. For the velocity criterion,
+the maximum value of ``dt_cfl`` is found and then a suitable timestep is found
+as::
+
+  dt_min_vel = hmin/max(dt_cfl)
+
+For the force based criterion we use the following::
+
+  dt_min_force = sqrt(hmin/sqrt(max(dt_force)))
+
+for the viscosity we have::
+
+  dt_min_visc = hmin/max(dt_visc_fac)
+
+Then the correct timestep is found as::
+
+  dt = cfl*min(dt_min_vel, dt_min_force, dt_min_visc)
+
+The ``cfl`` is set to 0.3 by default. One may pass ``--cfl`` to the
+application to change the CFL. Note that when the ``dt_adapt`` property is
+used the CFL has no effect as we assume that the user will compute a suitable
+value based on their requirements.
+
+The :py:class:`pysph.sph.integrator.Integrator` class code may be instructive
+to look at if you are wondering about any particular details.
 
 Illustration of the ``loop_all`` method
 ----------------------------------------

pysph/sph/integrator.py

@@ -43,6 +43,8 @@ def __init__(self, **kw):
         # This is set later when the underlying compiled integrator is created
         # by the SPHCompiler.
         self.c_integrator = None
+        self._has_dt_adapt = None
+        self.fixed_h = False
 
     def __repr__(self):
         name = self.__class__.__name__
@@ -71,6 +73,35 @@ def _get_dt_adapt_factors(self):
         cfl_f, force_f, visc_f = factors
         return cfl_f, force_f, visc_f
 
+    def _get_explicit_dt_adapt(self):
+        """Checks if the user is defining a 'dt_adapt' property where the
+        timestep is directly specified.
+
+        This returns None if no such parameter is found, else it returns the
+        allowed timestep.
+        """
+        a_eval = self.acceleration_evals[0]
+        if self._has_dt_adapt is None:
+            self._has_dt_adapt = any(
+                'dt_adapt' in pa.properties for pa in a_eval.particle_arrays
+            )
+        if self._has_dt_adapt:
+            dt_min = 1e20
+            for pa in a_eval.particle_arrays:
+                if 'dt_adapt' in pa.properties:
+                    if pa.gpu is not None:
+                        from compyle.array import minimum
+                        min_val = minimum(pa.gpu.dt_adapt)
+                    else:
+                        min_val = np.min(pa.dt_adapt)
+                    dt_min = min(dt_min, min_val)
+            if dt_min > 0.0:
+                return dt_min
+            else:
+                return None
+        else:
+            return None
+
     ##########################################################################
     # Public interface.
     ##########################################################################
@@ -117,6 +148,10 @@ def compute_time_step(self, dt, cfl):
         """If there are any adaptive timestep constraints, the appropriate
         timestep is returned, else None is returned.
         """
+        dt_adapt = self._get_explicit_dt_adapt()
+        if dt_adapt is not None:
+            return dt_adapt
+
         dt_cfl_fac, dt_force_fac, dt_visc_fac = self._get_dt_adapt_factors()
 
         # iterate over particles and find hmin if using variable h

pysph/sph/tests/test_integrator.py

@@ -105,6 +105,100 @@ def _integrate(self, integrator, dt, tf, post_step_callback):
             t += dt
 
 
+class EulerStep(IntegratorStep):
+    def stage1(self, d_idx, d_x, d_u, dt):
+        d_x[d_idx] += dt*d_u[d_idx]
+
+
+class TestIntegratorAdaptiveTimestep(TestIntegratorBase):
+    def test_compute_timestep_without_adaptive(self):
+        # Given.
+        integrator = EulerIntegrator(fluid=EulerStep())
+        equations = [SHM(dest="fluid", sources=None)]
+        self._setup_integrator(equations=equations, integrator=integrator)
+
+        # When
+        dt = integrator.compute_time_step(0.1, 0.5)
+
+        # Then
+        self.assertEqual(dt, None)
+
+    def test_compute_timestep_with_dt_adapt(self):
+        # Given.
+        self.pa.extend(1)
+        self.pa.align_particles()
+        self.pa.add_property('dt_adapt')
+        self.pa.dt_adapt[:] = [0.1, 0.2]
+
+        integrator = EulerIntegrator(fluid=EulerStep())
+        equations = [SHM(dest="fluid", sources=None)]
+        self._setup_integrator(equations=equations, integrator=integrator)
+
+        # When
+        dt = integrator.compute_time_step(0.1, 0.5)
+
+        # Then
+        self.assertEqual(dt, 0.1)
+
+    def test_compute_timestep_with_dt_adapt_with_invalid_values(self):
+        # Given.
+        self.pa.extend(1)
+        self.pa.align_particles()
+        self.pa.add_property('dt_adapt')
+        self.pa.dt_adapt[:] = [0.0, -2.0]
+
+        integrator = EulerIntegrator(fluid=EulerStep())
+        equations = [SHM(dest="fluid", sources=None)]
+        self._setup_integrator(equations=equations, integrator=integrator)
+
+        # When
+        dt = integrator.compute_time_step(0.1, 0.5)
+
+        # Then
+        self.assertEqual(dt, None)
+
+    def test_compute_timestep_with_dt_adapt_trumps_dt_cfl(self):
+        # Given.
+        self.pa.extend(1)
+        self.pa.align_particles()
+        self.pa.add_property('dt_adapt')
+        self.pa.add_property('dt_cfl')
+        self.pa.h[:] = 1.0
+        self.pa.dt_adapt[:] = [0.1, 0.2]
+        self.pa.dt_cfl[:] = 1.0
+
+        integrator = EulerIntegrator(fluid=EulerStep())
+        equations = [SHM(dest="fluid", sources=None)]
+        self._setup_integrator(equations=equations, integrator=integrator)
+
+        # When
+        cfl = 0.5
+        dt = integrator.compute_time_step(0.1, cfl)
+
+        # Then
+        self.assertEqual(dt, 0.1)
+
+    def test_compute_timestep_with_dt_cfl(self):
+        # Given.
+        self.pa.extend(1)
+        self.pa.align_particles()
+        self.pa.add_property('dt_cfl')
+        self.pa.h[:] = 1.0
+        self.pa.dt_cfl[:] = [1.0, 2.0]
+
+        integrator = EulerIntegrator(fluid=EulerStep())
+        equations = [SHM(dest="fluid", sources=None)]
+        self._setup_integrator(equations=equations, integrator=integrator)
+
+        # When
+        cfl = 0.5
+        dt = integrator.compute_time_step(0.1, cfl)
+
+        # Then
+        expect = cfl*1.0/(2.0)
+        self.assertEqual(dt, expect)
+
+
 class S1Step(IntegratorStep):
 
     def py_stage1(self, dest, t, dt):