Merged in prabhuramachandran/pysph/fix_tdamp (pull request #144)

Fix initial timestep damping and fix pre/post_stage callbacks

examples/SPHysics/dambreak_sphysics.py

@@ -57,7 +57,7 @@ def create_particles(ipart='IPART.gz', indat='INDAT.gz', **kwargs):
 
 # Create a solver. damping time is taken as 0.1% of the final time
 solver = Solver(dim=dim, kernel=kernel, integrator=integrator,
-                adaptive_timestep=True, tf=tf, dt=dt, tdamp=tf/1000)
+                adaptive_timestep=True, tf=tf, dt=dt, n_damp=100)
 
 # create the equations
 equations = [

examples/dam_break.py

@@ -130,9 +130,9 @@
 #integrator = EPECIntegrator(fluid=WCSPHStep(), boundary=WCSPHStep())
 integrator = TVDRK3Integrator(fluid=WCSPHTVDRK3Step(), boundary=WCSPHTVDRK3Step())
 
-# Create a solver. Damping time is taken as 0.1% of the final time
+# Create a solver.  The damping is performed for the first 50 iterations.
 solver = Solver(kernel=kernel, dim=dim, integrator=integrator,
-                dt=dt, tf=tf, adaptive_timestep=True, tdamp=tf/1000.0,
+                dt=dt, tf=tf, adaptive_timestep=True, n_damp=50,
                 fixed_h=False)
 
 # create the equations

examples/dam_break3D.py

@@ -120,7 +120,7 @@
 
 # Create a solver.
 solver = Solver(kernel=kernel, dim=dim, integrator=integrator, tf=tf, dt=dt,
-                adaptive_timestep=True, tdamp=tf/1000.0)
+                adaptive_timestep=True, n_damp=50)
 
 # create the equations
 equations = [

examples/elliptical_drop.py

@@ -115,18 +115,18 @@ def get_circular_patch(dx=0.025, **kwargs):
 # Create the Integrator. Currently, PySPH supports multi-stage,
 # predictor corrector and a TVD-RK3 integrators.
 
-#integrator = PECIntegrator(fluid=WCSPHStep()) 
+#integrator = PECIntegrator(fluid=WCSPHStep())
 #integrator = EPECIntegrator(fluid=WCSPHStep())
 integrator = TVDRK3Integrator( fluid=WCSPHTVDRK3Step() )
 
-# Construct the solver. tdamp determines the time until which smaller
-# time-steps are used when using adaptive time-steps. Use the tdamp
+# Construct the solver. n_damp determines the iterations until which smaller
+# time-steps are used when using adaptive time-steps. Use the output_at_times
 # list to specify instants of time at which the output solution is
 # required.
 dt = 5e-6; tf = 0.0075
 solver = Solver(kernel=kernel, dim=2, integrator=integrator,
                 dt=dt, tf=tf, adaptive_timestep=True,
-                cfl=0.05, tdamp=tf/1000.0,
+                cfl=0.05, n_damp=50,
                 output_at_times=[0.0033, 0.0052])
 
 # select True if you want to dump out remote particle properties in

examples/iisph/elliptical_drop.py

@@ -108,10 +108,8 @@ def get_circular_patch(dx=0.025, **kwargs):
 # Create the Integrator.
 integrator = EulerIntegrator(fluid=IISPHStep())
 
-# Construct the solver. tdamp determines the time until which smaller
-# time-steps are used when using adaptive time-steps. Use the output_at_times
-# list to specify instants of time at which the output solution is
-# required.
+# Construct the solver. Use the output_at_times list to specify instants of
+# time at which the output solution is  required.
 dt = 2e-4;
 tf = 0.0075
 solver = Solver(kernel=kernel, dim=2, integrator=integrator,

examples/transport_velocity/periodic_cylinders.py

@@ -143,10 +143,10 @@ def create_particles(**kwargs):
 # function evaluation per time step.
 integrator = PECIntegrator(fluid=TransportVelocityStep())
 
-# Create a solver. Damping time is taken as 0.1% of the final time
+# Create a solver. Damping is performed for 100 iterations.
 solver = Solver(
     kernel=kernel, dim=2, integrator=integrator,
-    adaptive_timestep=False, tf=tf, dt=dt, tdamp=tf/1000.0)
+    adaptive_timestep=False, tf=tf, dt=dt, n_damp=100)
 
 equations = [
 

pysph/solver/solver.py

@@ -9,7 +9,7 @@
 from pysph.sph.acceleration_eval import AccelerationEval
 from pysph.sph.sph_compiler import SPHCompiler
 
-from utils import FloatPBar, savez, load, dump
+from utils import FloatPBar, load, dump
 
 import logging
 logger = logging.getLogger(__name__)
@@ -29,9 +29,9 @@ class Solver(object):
 
     - t -- the internal time step counter
 
-    - pre_step_functions -- a list of functions to be performed before stepping
+    - pre_step_callbacks -- a list of functions to be performed before stepping
 
-    - post_step_functions -- a list of functions to execute after stepping
+    - post_step_callbacks -- a list of functions to execute after stepping
 
     - pfreq -- the output print frequency
 
@@ -44,7 +44,7 @@ class Solver(object):
     """
 
     def __init__(self, dim=2, integrator=None, kernel=None,
-                 tdamp=0.001, tf=1.0, dt=1e-3,
+                 n_damp=0, tf=1.0, dt=1e-3,
                  adaptive_timestep=False, cfl=0.3,
                  output_at_times = [],
                  fixed_h=False, **kwargs):
@@ -65,8 +65,9 @@ def __init__(self, dim=2, integrator=None, kernel=None,
         kernel : base.kernels.Kernel
             SPH kernel to use
 
-        tdamp : double
-            Time upto which damping of the initial solution is required
+        n_damp : int
+            Number of timesteps for which the initial damping is required.
+            Setting it to zero will disable damping the timesteps.
 
         tf, dt : double
             Final time and suggested initial time-step
@@ -77,7 +78,7 @@ def __init__(self, dim=2, integrator=None, kernel=None,
         cfl : double
             CFL number for adaptive time stepping
 
-        output_at_times : list
+        output_at_times : list/array
             Optional list of output times to force output
 
         fixed_h : bint
@@ -105,8 +106,8 @@ def __init__(self, dim=2, integrator=None, kernel=None,
         self.execute_commands = None
 
         # list of functions to be called before and after an integration step
-        self.pre_step_functions = []
-        self.post_step_functions = []
+        self.pre_step_callbacks = []
+        self.post_step_callbacks = []
 
         # List of functions to be called after each stage of the integrator.
         self.post_stage_callbacks = []
@@ -144,14 +145,14 @@ def __init__(self, dim=2, integrator=None, kernel=None,
         self.output_directory = self.fname+'_output'
 
         # solution damping to avoid impulsive starts
-        self.tdamp = tdamp
+        self.n_damp = n_damp
 
         # Use adaptive time steps and cfl number
         self.adaptive_timestep = adaptive_timestep
         self.cfl = cfl
 
         # list of output times
-        self.output_at_times = output_at_times
+        self.output_at_times = numpy.asarray(output_at_times)
         self.force_output = False
 
         # Use cell iterations or not.
@@ -168,6 +169,8 @@ def __init__(self, dim=2, integrator=None, kernel=None,
         # default time step constants
         self.tf = tf
         self.dt = dt
+        self._old_dt = dt
+        self._damping_factor = 1.0
 
         # flag for constant smoothing lengths
         self.fixed_h = fixed_h
@@ -226,6 +229,20 @@ def add_post_stage_callback(self, callback):
         """
         self.post_stage_callbacks.append(callback)
 
+    def add_post_step_callback(self, callback):
+        """These callbacks are called *after* each timestep is performed.
+
+        The callbacks are passed the solver instance (i.e. self).
+        """
+        self.post_step_callbacks.append(callback)
+
+    def add_pre_step_callback(self, callback):
+        """These callbacks are called *before* each timestep is performed.
+
+        The callbacks are passed the solver instance (i.e. self).
+        """
+        self.pre_step_callbacks.append(callback)
+
     def append_particle_arrrays(self, arrays):
         """ Append the particle arrays to the existing particle arrays
         """
@@ -309,7 +326,7 @@ def set_output_directory(self, path):
 
     def set_output_at_times(self, output_at_times):
         """ Set a list of output times """
-        self.output_at_times = output_at_times
+        self.output_at_times = numpy.asarray(output_at_times)
 
     def set_parallel_output_mode(self, mode="collected"):
         """Set the default solver dump mode in parallel.
@@ -353,9 +370,6 @@ def solve(self, show_progress=True):
         the stepping within the integrator.
 
         """
-        dt = self.dt
-        old_dt = self.dt
-
         if self.in_parallel:
             show = False
         else:
@@ -366,38 +380,36 @@ def solve(self, show_progress=True):
         self.dump_output()
         self.barrier() # everybody waits for this to complete
 
-        # initial solution damping time
-        tdamp = self.tdamp
-
         # Compute the accelerations once for the predictor corrector
         # integrator to work correctly at the first time step.
-        self.acceleration_eval.compute(self.t, dt)
+        self.acceleration_eval.compute(self.t, self.dt)
 
-        # solution output times
-        output_at_times = numpy.array( self.output_at_times )
+        # Now get a suitable adaptive (if requested) and damped timestep to
+        # integrate with.
+        self.dt = self._get_timestep()
 
         while self.t < self.tf:
 
             # perform any pre step functions
-            for func in self.pre_step_functions:
-                func.eval(self)
+            for callback in self.pre_step_callbacks:
+                callback(self)
 
             if self.rank == 0:
                 logger.debug(
                     "Iteration=%d, time=%f, timestep=%f" % \
-                        (self.count, self.t, dt)
+                        (self.count, self.t, self.dt)
                 )
             # perform the integration and update the time.
-            #print 'Solver Iteration', self.count, dt, self.t, tdamp
-            self.integrator.step(self.t, dt)
+            #print 'Solver Iteration', self.count, self.dt, self.t
+            self.integrator.step(self.t, self.dt)
 
             # perform any post step functions
-            for func in self.post_step_functions:
-                func.eval(self)
+            for callback in self.post_step_callbacks:
+                callback(self)
 
             # update time and iteration counters if successfully
             # integrated
-            self.t += dt
+            self.t += self.dt
             self.count += 1
 
             # dump output if the iteration number is a multiple of the
@@ -415,55 +427,16 @@ def solve(self, show_progress=True):
 
                 # re-set the time-step to the old time step before it
                 # was adjusted
-                dt = old_dt
+                self.dt = self._old_dt
 
             # update progress bar
             bar.update(self.t)
 
             # update the time for all arrays
             self.update_particle_time()
 
-            # compute the new time step across all processors
-            if self.adaptive_timestep:
-
-                # locally stable time step
-                dt = self.integrator.compute_time_step(
-                    self.dt, self.cfl)
-
-                # damp the initial solution
-                if self.t < tdamp:
-                    dt *= 0.5 * (numpy.sin(numpy.pi*(-0.5+self.t/tdamp)) + 1.0)
-
-                # set the globally stable time step
-                if self.in_parallel:
-                    dt = self.pm.update_time_steps(dt)
-
-            # adjust dt to land on final time
-            if self.t + dt > self.tf:
-                dt = self.tf - self.t
-
-            # adjust dt to land on specific output times
-            if output_at_times.size > 0:
-                tdiff = output_at_times - self.t
-                condition = (tdiff > 0) & (tdiff < dt)
-
-                if numpy.any( condition ):
-                    output_time = output_at_times[ numpy.where(condition) ]
-                    if abs(output_time - self.t) > 1e-14:
-                        # It sometimes happens that the current time is just
-                        # shy of the requested output time which results in a
-                        # ridiculously small dt so we skip that case.
-
-                        # save the old time-step and compute the new
-                        # time-step to fall on the specified output time
-                        # instant
-                        old_dt = dt
-                        dt = float( output_time - self.t )
-
-                        self.force_output = True
-
-            # set the adjusted time-step for the solver
-            self.dt = dt
+            # Compute the next timestep.
+            self.dt = self._get_timestep()
 
             if self.execute_commands is not None:
                 if self.count % self.command_interval == 0:
@@ -607,8 +580,81 @@ def setup_solver(self, options=None):
     ##########################################################################
     # Non-public interface.
     ##########################################################################
+    def _compute_timestep(self):
+        undamped_dt = self._get_undamped_timestep()
+        if self.adaptive_timestep:
+            # locally stable time step
+            dt = self.integrator.compute_time_step(undamped_dt, self.cfl)
+
+            # set the globally stable time step across all processors
+            if self.in_parallel:
+                dt = self.pm.update_time_steps(dt)
+        else:
+            dt = undamped_dt
+
+        return dt
+
+    def _damp_timestep(self, dt):
+        """Damp the timestep initially to prevent transient errors at startup.
+
+        This basically damps the initial timesteps by the factor
+
+        0.5  (sin(pi*(-0.5 + count/n_damp)) + 1)
+
+        Where n_damp is the number of iterations to damp the timestep for and
+        count is the number of iterations.
+
+        """
+        n_damp = self.n_damp
+        if self.count < n_damp and n_damp > 0:
+            iter_fraction = (self.count+1)/float(n_damp)
+            fac = 0.5*(numpy.sin(numpy.pi*(-0.5 + iter_fraction)) + 1.0)
+            self._damping_factor = fac
+        else:
+            self._damping_factor = 1.0
+
+        return dt*self._damping_factor
+
+    def _get_timestep(self):
+
+        dt = self._compute_timestep()
+        dt = self._damp_timestep(dt)
+
+        # adjust dt to land on final time
+        if self.t + dt > self.tf:
+            dt = self.tf - self.t
+
+        # solution output times
+        output_at_times = self.output_at_times
+
+        # adjust dt to land on specific output times
+        if len(output_at_times) > 0:
+            tdiff = output_at_times - self.t
+            condition = (tdiff > 0) & (tdiff < dt)
+
+            if numpy.any( condition ):
+                output_time = output_at_times[ numpy.where(condition) ]
+                if abs(output_time - self.t) > 1e-14:
+                    # It sometimes happens that the current time is just
+                    # shy of the requested output time which results in a
+                    # ridiculously small dt so we skip that case.
+
+                    # save the old time-step and compute the new
+                    # time-step to fall on the specified output time
+                    # instant
+                    self._old_dt = dt
+                    dt = float( output_time - self.t )
+
+                    self.force_output = True
+
+        return dt
+
+    def _get_undamped_timestep(self):
+        return self.dt/self._damping_factor
+
     def _post_stage_callback(self, time, dt, stage):
         for callback in self.post_stage_callbacks:
             callback(time, dt, stage)
 
+
 ############################################################################