[1D] Add callback hooks to make debugging 1D solver easier

This allows one to define functions in Python which will be called after each
successful time step or steady state solve, from which the state of the solver
can be inspected.

include/cantera/oneD/OneDim.h

@@ -301,6 +301,14 @@ public:
         m_interrupt = interrupt;
     }
 
+    //! Set a function that will be called after each successful timestep. The
+    //! function will be called with the size of the timestep as the argument.
+    //! Intended to be used for observing solver progress for debugging
+    //! purposes.
+    void setTimeStepCallback(Func1* callback) {
+        m_time_step_callback = callback;
+    }
+
 protected:
     void evalSSJacobian(doublereal* x, doublereal* xnew);
 
@@ -333,6 +341,9 @@ protected:
     //! Function called at the start of every call to #eval.
     Func1* m_interrupt;
 
+    //! User-supplied function called after each successful timestep.
+    Func1* m_time_step_callback;
+
     //! Number of time steps taken in the current call to solve()
     int m_nsteps;
 

include/cantera/oneD/Sim1D.h

@@ -169,6 +169,13 @@ public:
 
     virtual void resize();
 
+    //! Set a function that will be called after each successful steady-state
+    //! solve, before regridding. Intended to be used for observing solver
+    //! progress for debugging purposes.
+    void setSteadyCallback(Func1* callback) {
+        m_steady_callback = callback;
+    }
+
 protected:
     //! the solution vector
     vector_fp m_x;
@@ -194,6 +201,9 @@ protected:
     //! solution
     vector_int m_steps;
 
+    //! User-supplied function called after a successful steady-state solve.
+    Func1* m_steady_callback;
+
 private:
     /// Calls method _finalize in each domain.
     void finalize();

interfaces/cython/cantera/_cantera.pxd

@@ -720,6 +720,8 @@ cdef extern from "cantera/oneD/Sim1D.h":
         void setGridMin(int, double) except +
         void setFixedTemperature(double)
         void setInterrupt(CxxFunc1*) except +
+        void setTimeStepCallback(CxxFunc1*)
+        void setSteadyCallback(CxxFunc1*)
 
 cdef extern from "<sstream>":
     cdef cppclass CxxStringStream "std::stringstream":
@@ -1008,6 +1010,8 @@ cdef class Sim1D:
     cdef object _initial_guess_args
     cdef object _initial_guess_kwargs
     cdef Func1 interrupt
+    cdef Func1 time_step_callback
+    cdef Func1 steady_callback
 
 cdef class ReactionPathDiagram:
     cdef CxxReactionPathDiagram diagram

interfaces/cython/cantera/onedim.pyx

@@ -548,6 +548,28 @@ cdef class Sim1D:
         self.interrupt = f
         self.sim.setInterrupt(self.interrupt.func)
 
+    def set_time_step_callback(self, f):
+        """
+        Set a callback function to be called after each successful timestep.
+        The signature of *f* is `float f(float)`. The argument passed to *f* is
+        the size of the timestep. The output is ignored.
+        """
+        if not isinstance(f, Func1):
+            f = Func1(f)
+        self.time_step_callback = f
+        self.sim.setTimeStepCallback(self.time_step_callback.func)
+
+    def set_steady_callback(self, f):
+        """
+        Set a callback function to be called after each successful timestep.
+        The signature of *f* is `float f(float)`. The argument passed to *f* is
+        "0" and the output is ignored.
+        """
+        if not isinstance(f, Func1):
+            f = Func1(f)
+        self.steady_callback = f
+        self.sim.setSteadyCallback(self.steady_callback.func)
+
     def domain_index(self, dom):
         """
         Get the index of a domain, specified either by name or as a Domain1D

interfaces/cython/cantera/test/test_onedim.py

@@ -163,7 +163,6 @@ class TestFreeFlame(utilities.CanteraTest):
 
     def solve_mix(self, ratio=3.0, slope=0.3, curve=0.2, prune=0.0, refine=True):
         # Solve with the energy equation enabled
-
         self.sim.set_refine_criteria(ratio=ratio, slope=slope, curve=curve, prune=prune)
         self.sim.energy_enabled = True
         self.sim.solve(loglevel=0, refine_grid=refine)
@@ -530,9 +529,26 @@ class TestDiffusionFlame(utilities.CanteraTest):
         referenceFile = '../data/DiffusionFlameTest-h2-auto.csv'
         self.create_sim(p=ct.one_atm, mdot_fuel=2, mdot_ox=3)
 
+        nPoints = []
+        timesteps = []
+        def steady_func(x):
+            nPoints.append(len(self.sim.T))
+            return 0
+
+        def time_step_func(dt):
+            timesteps.append(dt)
+            self.assertGreater(dt, 0)
+            return 0
+
+        self.sim.set_steady_callback(steady_func)
+        self.sim.set_time_step_callback(time_step_func)
+
         self.sim.set_refine_criteria(ratio=3.0, slope=0.1, curve=0.12, prune=0.0)
         self.sim.solve(loglevel=0, auto=True)
 
+        self.assertNotEqual(len(nPoints), 0)
+        self.assertNotEqual(len(timesteps), 0)
+
         data = np.empty((self.sim.flame.n_points, self.gas.n_species + 4))
         data[:,0] = self.sim.grid
         data[:,1] = self.sim.u

src/oneD/OneDim.cpp

@@ -18,7 +18,8 @@ OneDim::OneDim()
       m_bw(0), m_size(0),
       m_init(false), m_pts(0), m_solve_time(0.0),
       m_ss_jac_age(20), m_ts_jac_age(20),
-      m_interrupt(0), m_nsteps(0), m_nsteps_max(500),
+      m_interrupt(0), m_time_step_callback(0),
+      m_nsteps(0), m_nsteps_max(500),
       m_nevals(0), m_evaltime(0.0)
 {
     m_newt.reset(new MultiNewton(1));
@@ -30,7 +31,8 @@ OneDim::OneDim(vector<Domain1D*> domains) :
     m_bw(0), m_size(0),
     m_init(false), m_solve_time(0.0),
     m_ss_jac_age(20), m_ts_jac_age(20),
-    m_interrupt(0), m_nsteps(0), m_nsteps_max(500),
+    m_interrupt(0), m_time_step_callback(0),
+    m_nsteps(0), m_nsteps_max(500),
     m_nevals(0), m_evaltime(0.0)
 {
     // create a Newton iterator, and add each domain.
@@ -223,6 +225,7 @@ int OneDim::solve(doublereal* x, doublereal* xnew, int loglevel)
         m_jac->updateTransient(m_rdt, m_mask.data());
         m_jac_ok = true;
     }
+
     return m_newt->solve(x, xnew, *this, *m_jac, loglevel);
 }
 
@@ -373,6 +376,9 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
             if (m == 100) {
                 dt *= 1.5;
             }
+            if (m_time_step_callback) {
+                m_time_step_callback->eval(dt);
+            }
             dt = std::min(dt, m_tmax);
             if (m_nsteps >= m_nsteps_max) {
                 throw CanteraError("OneDim::timeStep",

src/oneD/Sim1D.cpp

@@ -8,6 +8,7 @@
 #include "cantera/oneD/MultiNewton.h"
 #include "cantera/numerics/funcs.h"
 #include "cantera/base/xml.h"
+#include "cantera/numerics/Func1.h"
 
 using namespace std;
 
@@ -15,7 +16,8 @@ namespace Cantera
 {
 
 Sim1D::Sim1D(vector<Domain1D*>& domains) :
-    OneDim(domains)
+    OneDim(domains),
+    m_steady_callback(0)
 {
     // resize the internal solution vector and the work array, and perform
     // domain-specific initialization of the solution vector.
@@ -257,6 +259,10 @@ void Sim1D::solve(int loglevel, bool refine_grid)
                     }
                     writelog("] point grid(s).\n");
                 }
+                if (m_steady_callback) {
+                    m_steady_callback->eval(0);
+                }
+
                 if (loglevel > 6) {
                     save("debug_sim1d.xml", "debug",
                          "After successful Newton solve");