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");