Fixed #56: Can specify no. of tolerated errors in solver.run.

cmf/__init__.py

@@ -22,8 +22,8 @@
 from .stopwatch import StopWatch
 
 
-__version__ = '1.4.2a'
-__compiletime__ = 'Tue Mar 12 21:43:31 2019'
+__version__ = '1.5.0a'
+__compiletime__ = 'Thu Mar 14 23:07:32 2019'
 
 from .cmf_core import connect_cells_with_flux as __ccwf
 

cmf/cmf_core.py

@@ -13397,22 +13397,88 @@ def __len__(self, *args, **kwargs):
         return _cmf_core.Integrator___len__(self, *args, **kwargs)
 
 
-    def __call__(self,t,reset=False):
-        if t<self.t:
-            self.integrate_until(self.t+t,reset=reset)
-        else:
-            self.integrate_until(t,reset=reset)
-        return self.t
     t = property(get_t,set_t,doc="Sets the actual time of the solution")
     dt = property(get_dt,doc="Get the current time step of the solver")
-    def run(self,start=None,end=None,step=day*1):
+    def __call__(self, t, dt=None, reset=False):
+        """
+        Advances the integration until `t`
+
+        A shortcut to .integrate_until
+
+        Parameters
+        ----------
+        t : cmf.Time
+            The time step to advance to. If t < current time, the solver will
+            advance to self.t + t
+        dt : cmf.Time, optional
+            The timestep for the integration. If not given try to integrate in one step
+        reset : bool, optional
+            If True, the solver will perform a reset before starting
+
+        Returns
+        -------
+        cmf.Time
+            The new time stamp
+        """
+        if dt is None:
+            dt = Time()
+        if t < self.t:
+            self.integrate_until(self.t+t, dt, reset=reset)
+        else:
+            self.integrate_until(t, dt, reset=reset)
+        return self.t
+
+    def run(self, start=None, end=None, step=day*1, max_errors=0, reset=False):
+        """
+        Returns an iterator over the timesteps start..end
+
+        **Examples:**
+
+        >>> solver=cmf.CVodeIntegrator(...)
+        >>> for t in solver.run(solver.t, solver.t + cmf.week, cmf.h):
+        >>>    print(t, solver[0].state)
+        or with list comprehension
+        >>> states = [solver[0].state for t in solver.run(solver.t, solver.t + cmf.week, cmf.h)]
+
+        Parameters
+        ----------
+        start : cmf.Time, optional
+                Start time for the solver iteration
+        end : cmf.Time, optional
+                End time of the iteration
+        step : cmf.Time, optional
+                Step size for the integration
+        max_errors: int
+                Number of tolerated errors. If >0, up to these number of runtime errors
+                will be saved with their time and the integration proceeds after a reset
+                of the solver. Some systems operate with values close to their physical
+                limits and inifinite values in the integration can easily occur. For
+                these kind of systems set max_errors to eg. 10. A larger number of errors
+                should be eliminated usually.
+        reset: bool
+                If True, the solver performs a `reset` at every time step
+
+        Yields
+        ------
+        cmf.Time
+             the actual timestep
+        """
+        import logging
         if not start is None:
-            self.t=start
+            self.t = start
         if end is None:
             end = self.t + 100*step
-        while self.t<end:
-            self(self.t+step)
-            yield self.t
+        errors = []
+        t = self.t
+        while self.t < end:
+            try:
+                t = self(self.t+step, reset=reset)
+            except Exception as e:
+                if len(errors) < max_errors:
+                    errors.append((t, e))
+                    self.reset()
+                    logging.warning(str(t) + ': ' + str(e))
+            yield t
 
 Integrator.__getitem__ = new_instancemethod(_cmf_core.Integrator___getitem__, None, Integrator)
 Integrator.get_dxdt = new_instancemethod(_cmf_core.Integrator_get_dxdt, None, Integrator)

cmf/cmf_core_src/math/integrator.i

@@ -31,21 +31,87 @@
         return $self->size();
     }
 %pythoncode {
-    def __call__(self,t,reset=False):
-        if t<self.t:
-            self.integrate_until(self.t+t,reset=reset)
-        else:
-            self.integrate_until(t,reset=reset)
-        return self.t
     t = property(get_t,set_t,doc="Sets the actual time of the solution")
     dt = property(get_dt,doc="Get the current time step of the solver")
-    def run(self,start=None,end=None,step=day*1):
+    def __call__(self, t, dt=None, reset=False):
+        """
+        Advances the integration until `t`
+
+        A shortcut to .integrate_until
+
+        Parameters
+        ----------
+        t : cmf.Time
+            The time step to advance to. If t < current time, the solver will
+            advance to self.t + t
+        dt : cmf.Time, optional
+            The timestep for the integration. If not given try to integrate in one step
+        reset : bool, optional
+            If True, the solver will perform a reset before starting
+
+        Returns
+        -------
+        cmf.Time
+            The new time stamp
+        """
+        if dt is None:
+            dt = Time()
+        if t < self.t:
+            self.integrate_until(self.t+t, dt, reset=reset)
+        else:
+            self.integrate_until(t, dt, reset=reset)
+        return self.t
+
+    def run(self, start=None, end=None, step=day*1, max_errors=0, reset=False):
+        """
+        Returns an iterator over the timesteps start..end
+
+        **Examples:**
+
+        >>> solver=cmf.CVodeIntegrator(...)
+        >>> for t in solver.run(solver.t, solver.t + cmf.week, cmf.h):
+        >>>    print(t, solver[0].state)
+        or with list comprehension
+        >>> states = [solver[0].state for t in solver.run(solver.t, solver.t + cmf.week, cmf.h)]
+
+        Parameters
+        ----------
+        start : cmf.Time, optional
+                Start time for the solver iteration
+        end : cmf.Time, optional
+                End time of the iteration
+        step : cmf.Time, optional
+                Step size for the integration
+        max_errors: int
+                Number of tolerated errors. If >0, up to these number of runtime errors
+                will be saved with their time and the integration proceeds after a reset
+                of the solver. Some systems operate with values close to their physical
+                limits and inifinite values in the integration can easily occur. For
+                these kind of systems set max_errors to eg. 10. A larger number of errors
+                should be eliminated usually.
+        reset: bool
+                If True, the solver performs a `reset` at every time step
+
+        Yields
+        ------
+        cmf.Time
+             the actual timestep
+        """
+        import logging
         if not start is None:
-            self.t=start
+            self.t = start
         if end is None:
             end = self.t + 100*step
-        while self.t<end:
-            self(self.t+step)
-            yield self.t
+        errors = []
+        t = self.t
+        while self.t < end:
+            try:
+                t = self(self.t+step, reset=reset)
+            except Exception as e:
+                if len(errors) < max_errors:
+                    errors.append((t, e))
+                    self.reset()
+                    logging.warning(str(t) + ': ' + str(e))
+            yield t
 }
 }

setup.py

@@ -30,7 +30,7 @@
 from distutils.sysconfig import customize_compiler
 from distutils.command.build_py import build_py
 
-version = '1.4.2a'
+version = '1.5.0a'
 branchversion = version
 try:
     from pygit2 import Repository

test/solvers.py

@@ -3,7 +3,7 @@
 
 """
 import cmf
-import numpy as np
+
 import unittest
 
 
@@ -79,11 +79,11 @@ def test_solver_results_solute(self):
             vol = [s.volume for s in stores]
             smass = [s[X].state for s in stores]
 
-            rmse_v = sum((vr - v)**2 for v, vr in zip(vol, vol_ref)) / len(vol)
-            rmse_s = sum((sr - s)**2 for s, sr in zip(smass, smass_ref)) / len(smass)
+            mse_v = sum((vr - v)**2 for v, vr in zip(vol, vol_ref)) / len(vol)
+            mse_s = sum((sr - s)**2 for s, sr in zip(smass, smass_ref)) / len(smass)
 
-            self.assertAlmostEqual(1 - rmse_v, 1, 2, "RMSE between reference volume and {} too large".format(st.__name__))
-            self.assertAlmostEqual(1 - rmse_s, 1, 2, "RMSE between reference solute and {} too large".format(st.__name__))
+            self.assertAlmostEqual(1 - mse_v, 1, 2, "RMSE between reference volume and {} too large".format(st.__name__))
+            self.assertAlmostEqual(1 - mse_s, 1, 2, "RMSE between reference solute and {} too large".format(st.__name__))
 
     def test_solver_results_nosolute(self):
 
@@ -95,14 +95,33 @@ def test_solver_results_nosolute(self):
             p, stores, X = get_project(False)
 
             solver = st(p)
-            solver.integrate_until(cmf.day * 3, cmf.h)
+            solver(cmf.day * 3, cmf.h)
 
             vol = [s.volume for s in stores]
 
-            rmse_v = sum((vr - v)**2 for v, vr in zip(vol, vol_ref)) / len(vol)
+            mse_v = sum((vr - v)**2 for v, vr in zip(vol, vol_ref)) / len(vol)
+
+            self.assertAlmostEqual(1 - mse_v, 1, 2, "MSE between reference volume and {} too large".format(st.__name__))
+
+    def test_solver_run(self):
+
 
-            self.assertAlmostEqual(1 - rmse_v, 1, 2, "RMSE between reference volume and {} too large".format(st.__name__))
+        for st in solver_types:
+
+            p, stores, X = get_project(False)
 
+            solver = st(p)
+            # Test all run parameters
+            for t in solver.run(cmf.Time(), cmf.day, cmf.h, reset=False, max_errors=2):
+                ...
+            self.assertEqual(solver.t, cmf.day)
+            # Test little parameter count
+            for t in solver.run(step=cmf.h):
+                ...
+            self.assertEqual(solver.t, cmf.day + 100 * cmf.h)
 # Check solver results
 # Check jacobian (CVodeDirect, CVodeKLU)
 
+
+if __name__ == '__main__':
+    unittest.main(verbosity=5)

tools/Doxyfile

@@ -38,7 +38,7 @@ PROJECT_NAME           = cmf
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = 1.4.2a
+PROJECT_NUMBER         = 1.5.0a
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a