dump times implementation

examples/example_utils.py

@@ -521,8 +521,8 @@ class HysopArgParser(argparse.ArgumentParser):
                 dest='dump_freq',
                 help=('HDF5 output frequency in terms of iterations (default=10).' 
                      +' Use 0 to disable frequency based dumping.'))
-        file_io.add_argument('--dump-times', type=str, default=None, 
-                action=self.split, container=tuple, append=True,
+        file_io.add_argument('--dump-times', type=str, default=None, convert=float,
+                action=self.split, container=tuple, append=False,
                 dest='dump_times',
                 help='Comma delimited list of additional HDF5 output times of interest.')
         file_io.add_argument('--cache-dir', type=str, default=None, 
@@ -552,13 +552,19 @@ class HysopArgParser(argparse.ArgumentParser):
             msg=msg.format(args.tstart, args.tend)
             self.error(msg)
 
+        args.dump_times = set(args.dump_times)
+
         msg='args.dump_times = {}\n'.format(args.dump_times)
         for dp in args.dump_times:
             if (dp<args.tstart):
                 msg+='Dump time of interest t={} happens before tstart={}.'
                 msg=msg.format(dp, args.tstart)
                 self.error(msg)
-            elif (dp>args.tend):
+            elif (dp==args.tend):
+                msg+='Dump time of interest t={} happens exactly at tend={}.'
+                msg=msg.format(dp, args.tend)
+                self.error(msg)
+            elif (dp>=args.tend):
                 msg+='Dump time of interest t={} happens after tend={}.'
                 msg=msg.format(dp, args.tend)
                 self.error(msg)

examples/shear_layer/shear_layer.py

@@ -147,8 +147,8 @@ def compute(args):
     simu = Simulation(start=args.tstart, end=args.tend, 
                       nb_iter=args.nb_iter,
                       max_iter=args.max_iter,
-                      dt0=args.dt0, 
-                        t=t, dt=dt)
+                      dt0=args.dt0, times_of_interest=args.dump_times,
+                      t=t, dt=dt)
     simu.write_parameters(t, dt, filename='parameters.txt', precision=4)
     
     # Initialize only the vorticity
@@ -250,9 +250,11 @@ if __name__=='__main__':
 
     parser.set_defaults(impl='cl', ndim=2, npts=(257,),
                         box_origin=(0.0,), box_length=(1.0,), 
-                        tstart=0.0, tend=1.20, 
+                        tstart=0.0, tend=1.25, 
                         nb_iter=None, dt=None, 
                         dt0=1e-4, cfl=0.5, lcfl=0.125,
-                        case=1, dump_times=(0.8, 1.20))
+                        case=1, 
+                        dump_freq=0,
+                        dump_times=(0.8, 1.20))
 
     parser.run(compute)

hysop/operator/hdf_io.py

@@ -258,14 +258,12 @@ class HDF_Writer(HDF_IO):
     def apply(self, simulation=None, **kwds):
         if simulation is None:
             raise ValueError("Missing simulation value for monitoring.")
-        if (self.io_params.frequency>0): 
-            ite = simulation.current_iteration
-            should_dump = (ite == -1)
-            should_dump |= (ite % self.io_params.frequency == 0)
-            should_dump |= simulation._next_is_last
-            if should_dump:
-                self.step(simulation)
-                self._count += 1
+        ite = simulation.current_iteration
+        should_dump  = (self.io_params.frequency>0) and (ite % self.io_params.frequency == 0)
+        should_dump |= simulation.is_time_of_interest
+        if should_dump:
+            self.step(simulation)
+            self._count += 1
 
     def _setup_grid_template(self):
         topo = self.topology
@@ -360,14 +358,14 @@ class HDF_Writer(HDF_IO):
         
         # Collect datas required to write the xdmf file
         # --> add tuples (counter, time).
-        if (simu.time == self._last_written_time):
+        if (simu.t() == self._last_written_time):
             msg = 'You cannot write two hdf files for the same '
             msg += '(time, var) set. '
             msg += 'If you want to save a field two times for '
             msg += 'a single time value, please use two hdf_writer operators.'
             raise RuntimeError(msg)
-        self._xdmf_data_files.append((self._count, simu.time))
-        self._last_written_time = simu.time
+        self._xdmf_data_files.append((self._count, simu.t()))
+        self._last_written_time = simu.t()
 
         self._hdf_file.close()
 

hysop/simulation.py

@@ -33,7 +33,7 @@ from hysop import dprint, vprint
 from hysop.deps import sys, os
 from hysop.constants import HYSOP_REAL
 from hysop.parameters.scalar_parameter import ScalarParameter
-from hysop.tools.types import first_not_None
+from hysop.tools.types import first_not_None, to_set
 from hysop.tools.numpywrappers import npw
 from hysop.tools.io_utils import IO, IOParams
 
@@ -46,7 +46,8 @@ class Simulation(object):
     """
 
     def __init__(self, name=None, start=0.0, end=1.0, nb_iter=None, dt0=None,
-                 max_iter=None, t=None, dt=None, **kwds):
+                 max_iter=None, t=None, dt=None, times_of_interest=None,
+                 **kwds):
         """
         Parameters
         ----------
@@ -67,6 +68,13 @@ class Simulation(object):
         max_iter : int, optional
             Maximum number of iterations allowed.
             Defaults to 1e9.
+        times_of_interest: array-like of float
+            List of times ti where the simulation may
+            modify current timestep to get t=ti.
+            Mainly used by HDF_Writers for precise 
+            time dependent dumping.
+            tstart < ti <= tend
+            Defaults to empty set.
 
         Attributes
         ----------
@@ -82,6 +90,8 @@ class Simulation(object):
             The scalar parameter that represents time.
         time_step: double
             Value of the dt parameter.
+        time_of_interest: float
+            Current simulation time target.
 
         Notes
         -----
@@ -96,17 +106,12 @@ class Simulation(object):
           with self.tkp1 == self.time
 
         """
-        # Simulation final time
         self.end = end
-        # Starting time
         self.start = start
-        # Simulation current time
         self.time = start
-        # Is simulation is terminated
         self.is_over = False
-        # Iteration counter
         self.current_iteration = -1
-        # Number of iterations
+
         if (nb_iter is not None):
             self.nb_iter = nb_iter
             msg = '------------------------------------------------\n'
@@ -150,19 +155,29 @@ class Simulation(object):
             assert isinstance(t, ScalarParameter), type(t)
             assert not t.const, 't cannot be a constant parameter.'
             t.value = start
-        # t + dt
+        
+        # tk+1 = t + dt
         self.tkp1 = start + self.time_step
-        assert self.end > self.start, \
+
+        assert (self.end > self.start), \
             'Final time must be greater than initial time'
         assert (self.start + self.time_step) <= self.end,\
             'start + step is bigger than end.'
 
+        # times of interest
+        times_of_interest = to_set(first_not_None(times_of_interest, []))
+        times_of_interest = tuple(sorted(times_of_interest))
+        for toi in times_of_interest:
+            assert self.start <= toi < self.end, toi
+        self.times_of_interest = times_of_interest
+
         # Internal tolerance for timer
         self.tol = eps
         # True if initialize has been called.
-        self._is_ready = False
-        self._next_is_last = False
-        self._parameters_to_write = []
+        self._is_ready                 = False
+        self._next_is_last             = False
+        self._next_is_time_of_interest = False
+        self._parameters_to_write      = []
 
     def _get_time_step(self):
         """Get current timestep."""
@@ -201,26 +216,57 @@ class Simulation(object):
             return
 
         self.t.set_value(self.tkp1)
+        
+        self.is_time_of_interest = False
+        if (self.target_time_of_interest is not None):
+            if (abs(self.tkp1 - self.target_time_of_interest) <= self.tol):
+                self.next_time_of_interest()
+                self.is_time_of_interest = True
+        
         self.tkp1 = self.t() + self.time_step
         if abs(self.tkp1 - self.end) <= self.tol:
             self._next_is_last = True
-        elif self.tkp1 > self.end:
-            # resize ...
-            dprint('> Next iteration is last iteration, clamping dt.')
+        elif (self.tkp1 > self.end):
+            msg='** Next iteration is last iteration, clamping dt to achieve t={}. **'
+            msg=msg.format(self.end)
+            vprint()
+            self._print_banner(msg)
             self._next_is_last = True
             self.tkp1 = self.end
             self.update_time_step(self.end - self.t())
+        elif (self.target_time_of_interest is not None) and \
+                (self.tkp1 > self.target_time_of_interest):
+            msg='** Next iteration is a time of interest, clamping dt to achieve t={}. **'
+            msg=msg.format(self.target_time_of_interest)
+            vprint()
+            self._print_banner(msg)
+            self.tkp1 = self.target_time_of_interest
+            self.update_time_step(self.target_time_of_interest - self.t())
+
         self.current_iteration += 1
         self.time = self.tkp1
 
-        if self.current_iteration + 2 > self.max_iter:
+        if (self.current_iteration + 2 > self.max_iter):
             msg = '** Next iteration will be the last because max_iter={} will be achieved. **'
             msg=msg.format(self.max_iter)
             vprint()
-            vprint('*'*len(msg))
-            vprint(msg)
-            vprint('*'*len(msg))
+            self._print_banner(msg)
             self._next_is_last = True
+            self.is_time_of_interest = True
+
+    def _print_banner(self, msg):
+        vprint('*'*len(msg))
+        vprint(msg)
+        vprint('*'*len(msg))
+
+    def next_time_of_interest(self):
+        toi_counter       = self.toi_counter
+        times_of_interest = self.times_of_interest
+        if (toi_counter<len(times_of_interest)):
+            self.target_time_of_interest = times_of_interest[toi_counter]
+            self.toi_counter += 1
+        else:
+            self.target_time_of_interest = None
 
     def update_time_step(self, dt):
         """Update time step for the next iteration
@@ -243,11 +289,27 @@ class Simulation(object):
         """(Re)set simulation to initial values
         --> back to iteration 0 and ready to run.
         """
-        self.t.set_value(self.start)
-        self.update_time_step(self._dt0)
-        self.tkp1 = self.start + self.time_step
-        assert self.tkp1 <= self.end
-        if abs(self._dt0 - self.end) <= self.tol:
+        tstart, tend = self.start, self.end
+        times_of_interest = self.times_of_interest
+
+        self.toi_counter = 0
+        self.next_time_of_interest()
+        self.is_time_of_interest = False
+
+        assert ((tend - tstart) >= self.tol)
+        assert (tend >= self.target_time_of_interest >= tstart)
+        if abs(self.target_time_of_interest - tstart) <= self.tol:
+            self.next_time_of_interest()
+            self.is_time_of_interest = True
+        assert (tend >= self.target_time_of_interest > tstart)
+
+        dt0 = min(self._dt0, self.target_time_of_interest-tstart)
+        self.t.set_value(tstart)
+        self.update_time_step(dt0)
+        self.tkp1 = tstart + self.time_step
+        assert self.tkp1 < self.target_time_of_interest <= tend
+
+        if abs(self.tkp1 - self.end) <= self.tol:
             self._next_is_last = True
         else:
             self._next_is_last = False