From efe43b389e2a79eaa45434e99c432c0ce690abf6 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 14 May 2016 17:24:16 -0400 Subject: [PATCH] [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 | 11 ++++++++++ include/cantera/oneD/Sim1D.h | 10 +++++++++ interfaces/cython/cantera/_cantera.pxd | 4 ++++ interfaces/cython/cantera/onedim.pyx | 22 +++++++++++++++++++ interfaces/cython/cantera/test/test_onedim.py | 18 ++++++++++++++- src/oneD/OneDim.cpp | 10 +++++++-- src/oneD/Sim1D.cpp | 8 ++++++- 7 files changed, 79 insertions(+), 4 deletions(-) diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index 15db13260..8b32af890 100644 --- a/include/cantera/oneD/OneDim.h +++ b/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; diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index 928624a76..a5bb5b05d 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/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(); diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 7e04cae5f..54f5198bf 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/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 "": 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 diff --git a/interfaces/cython/cantera/onedim.pyx b/interfaces/cython/cantera/onedim.pyx index 6fd7729bc..b88dc6afe 100644 --- a/interfaces/cython/cantera/onedim.pyx +++ b/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 diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py index 009049a84..4f46bc741 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/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 diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index dc75404e8..e3203fa98 100644 --- a/src/oneD/OneDim.cpp +++ b/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 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", diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index 8eaae6fcd..8b0748a8b 100644 --- a/src/oneD/Sim1D.cpp +++ b/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& 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");