From 848a3bf0e329034d5c5d5ba1b2d8c846ac725ef4 Mon Sep 17 00:00:00 2001 From: Nick Date: Thu, 7 Mar 2019 11:05:39 -0500 Subject: [PATCH] Add InterfaceKinetics.advance_coverages integrator options to the cython interface, and test --- include/cantera/kinetics/ImplicitSurfChem.h | 34 ++++++++++-- include/cantera/kinetics/InterfaceKinetics.h | 14 ++--- interfaces/cython/cantera/_cantera.pxd | 2 +- interfaces/cython/cantera/kinetics.pyx | 9 ++-- interfaces/cython/cantera/test/test_onedim.py | 2 +- .../cython/cantera/test/test_reactor.py | 37 +++++++++++++ src/kinetics/ImplicitSurfChem.cpp | 53 +++++++++++++------ src/kinetics/InterfaceKinetics.cpp | 8 +-- 8 files changed, 124 insertions(+), 35 deletions(-) diff --git a/include/cantera/kinetics/ImplicitSurfChem.h b/include/cantera/kinetics/ImplicitSurfChem.h index 1ab6f7c0b..3ad6c774b 100644 --- a/include/cantera/kinetics/ImplicitSurfChem.h +++ b/include/cantera/kinetics/ImplicitSurfChem.h @@ -63,14 +63,17 @@ public: * freedom representing the concentration of surface adsorbates. * @param rtol The relative tolerance for the integrator * @param atol The absolute tolerance for the integrator - * @param maxStepSize The maximum step-size the integrator is allowed to take - * @param maxSteps the maximum number of time-steps the integrator can take + * @param maxStepSize The maximum step-size the integrator is allowed to take. + * If zero, this option is disabled. + * @param maxSteps The maximum number of time-steps the integrator can take. + * If not supplied, uses the default value in the CVodesIntegrator (20000). * @param maxErrTestFails the maximum permissible number of error test failures + * If not supplied, uses the default value in CVODES (7). */ ImplicitSurfChem(std::vector k, - doublereal rtol=1.e-7, doublereal atol=1.e-14, - doublereal maxStepSize=0, size_t maxSteps=0, - size_t maxErrTestFails=0); + double rtol=1.e-7, double atol=1.e-14, + double maxStepSize=0, size_t maxSteps=20000, + size_t maxErrTestFails=7); virtual ~ImplicitSurfChem() {}; @@ -79,6 +82,27 @@ public: */ virtual void initialize(doublereal t0 = 0.0); + /*! + * Set the maximum integration step-size. Note, setting this value to zero + * disables this option + */ + virtual void setMaxStepSize(double maxstep = 0.0); + + /*! + * Set the relative and absolute integration tolerances. + */ + virtual void setTolerances(double rtol=1.e-7, double atol=1.e-14); + + /*! + * Set the maximum number of CVODES integration steps. + */ + virtual void setMaxSteps(size_t maxsteps = 20000); + + /*! + * Set the maximum number of CVODES error test failures + */ + virtual void setMaxErrTestFails(size_t maxErrTestFails = 7); + //! Integrate from t0 to t1. The integrator is reinitialized first. /*! * This routine does a time accurate solve from t = t0 to t = t1. diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index bf9823f05..0b93fb4e8 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -242,14 +242,16 @@ public: * @param tstep Time value to advance the surface coverages * @param rtol The relative tolerance for the integrator * @param atol The absolute tolerance for the integrator - * @param maxStepSize The maximum step-size the integrator is allowed to take - * @param maxSteps the maximum number of time-steps the integrator can take - * before reaching tstep + * @param maxStepSize The maximum step-size the integrator is allowed to take. + * If zero, this option is disabled. + * @param maxSteps The maximum number of time-steps the integrator can take. + * If not supplied, uses the default value in CVodeIntegrator (20000). * @param maxErrTestFails the maximum permissible number of error test failures + * If not supplied, uses the default value in CVODES (7). */ - void advanceCoverages(doublereal tstep, doublereal rtol=1.e-7, - doublereal atol=1.e-14, doublereal maxStepSize=0, - size_t maxSteps=0, size_t maxErrTestFails=0); + void advanceCoverages(doublereal tstep, double rtol=1.e-7, + double atol=1.e-14, double maxStepSize=0, + size_t maxSteps=20000, size_t maxErrTestFails=7); //! Solve for the pseudo steady-state of the surface problem /*! diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 6db718c55..03992e958 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -380,7 +380,7 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera": cdef extern from "cantera/kinetics/InterfaceKinetics.h": cdef cppclass CxxInterfaceKinetics "Cantera::InterfaceKinetics": - void advanceCoverages(double) except +translate_exception + void advanceCoverages(double, double, double, double, size_t, size_t) except +translate_exception void solvePseudoSteadyStateProblem() except +translate_exception diff --git a/interfaces/cython/cantera/kinetics.pyx b/interfaces/cython/cantera/kinetics.pyx index 1eed3a12b..0c7ee9170 100644 --- a/interfaces/cython/cantera/kinetics.pyx +++ b/interfaces/cython/cantera/kinetics.pyx @@ -376,13 +376,16 @@ cdef class InterfaceKinetics(Kinetics): self._phase_indices[phase.name] = i self._phase_indices[i] = i - def advance_coverages(self, double dt): + def advance_coverages(self, double dt, double rtol=1e-7, double atol=1e-14, + double max_step_size=0, int max_steps=20000, + int max_error_test_failures=7): """ This method carries out a time-accurate advancement of the surface coverages for a specified amount of time. """ - (self.kinetics).advanceCoverages(dt) - + (self.kinetics).advanceCoverages( + dt, rtol, atol, max_step_size, max_steps, max_error_test_failures) + def advance_coverages_to_steady_state(self): """ This method advances the surface coverages to steady state. diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py index 2690f0c59..778fa8e3b 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/interfaces/cython/cantera/test/test_onedim.py @@ -896,7 +896,7 @@ class TestImpingingJet(utilities.CanteraTest): # integrate the coverage equations holding the gas composition fixed # to generate a good starting estimate for the coverages. - surf_phase.advance_coverages(1.0) + surf_phase.advance_coverages(1.) sim = ct.ImpingingJet(gas=gas, width=width, surface=surf_phase) sim.set_refine_criteria(10.0, 0.3, 0.4, 0.0) diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py index 11f759a17..2e0e67f09 100644 --- a/interfaces/cython/cantera/test/test_reactor.py +++ b/interfaces/cython/cantera/test/test_reactor.py @@ -1542,3 +1542,40 @@ class PureFluidReactorTest(utilities.CanteraTest): self.assertEqual(states.X[-2], 1) for i in range(3,7): self.assertNear(states.T[i], states.T[2]) + + +class AdvanceCoveragesTest(utilities.CanteraTest): + def setup(self, model='ptcombust.xml', gas_phase='gas', + interface_phase='Pt_surf'): + # create gas and interface + self.gas = ct.Solution('ptcombust.xml', 'gas') + self.surf = ct.Interface('ptcombust.xml', 'Pt_surf', [self.gas]) + + def test_advance_coverages_parameters(self): + # create gas and interface + self.setup() + + # first, test max step size & max steps + dt = 1.0 + max_steps = 10 + max_step_size = dt / (max_steps + 1) + # this should throw an error, as we can't reach dt + with self.assertRaises(ct.CanteraError): + self.surf.advance_coverages( + dt=dt, max_step_size=max_step_size, max_steps=max_steps) + + # next, run with different tolerances + self.setup() + self.surf.coverages = 'O(S):0.1, PT(S):0.5, H(S):0.4' + self.gas.TP = self.surf.TP + + self.surf.advance_coverages(dt=dt, rtol=1e-5, atol=1e-12) + cov = self.surf.coverages[:] + + self.surf.coverages = 'O(S):0.1, PT(S):0.5, H(S):0.4' + self.gas.TP = self.surf.TP + self.surf.advance_coverages(dt=dt, rtol=1e-7, atol=1e-14) + + # check that the solutions are similar, but not identical + self.assertArrayNear(cov, self.surf.coverages) + self.assertTrue(any(cov != self.surf.coverages)) diff --git a/src/kinetics/ImplicitSurfChem.cpp b/src/kinetics/ImplicitSurfChem.cpp index 1586e72b2..771a34330 100644 --- a/src/kinetics/ImplicitSurfChem.cpp +++ b/src/kinetics/ImplicitSurfChem.cpp @@ -18,8 +18,8 @@ namespace Cantera { ImplicitSurfChem::ImplicitSurfChem( - vector k, doublereal rtol, doublereal atol, - doublereal maxStepSize, size_t maxSteps, + vector k, double rtol, double atol, + double maxStepSize, size_t maxSteps, size_t maxErrTestFails) : m_nv(0), m_numTotalBulkSpecies(0), @@ -109,28 +109,49 @@ void ImplicitSurfChem::getState(doublereal* c) } } -void ImplicitSurfChem::initialize(doublereal t0) +void ImplicitSurfChem::setMaxStepSize(double maxstep) { - m_integ->setTolerances(m_rtol, m_atol); - if (m_maxstep > 0) - { + m_maxstep = maxstep; + if (m_maxstep > 0) { m_integ->setMaxStepSize(m_maxstep); } - if (m_nmax > 0) - { - m_integ->setMaxSteps(m_nmax); - } - if (m_maxErrTestFails > 0) - { - m_integ->setMaxErrTestFails(m_maxErrTestFails); - } +} + +void ImplicitSurfChem::setTolerances(double rtol, double atol) +{ + m_rtol = rtol; + m_atol = atol; + m_integ->setTolerances(m_rtol, m_atol); +} + +void ImplicitSurfChem::setMaxSteps(size_t maxsteps) +{ + m_nmax = maxsteps; + m_integ->setMaxSteps(m_nmax); +} + +void ImplicitSurfChem::setMaxErrTestFails(size_t maxErrTestFails) +{ + m_maxErrTestFails = maxErrTestFails; + m_integ->setMaxErrTestFails(m_maxErrTestFails); +} + +void ImplicitSurfChem::initialize(doublereal t0) +{ + this->setTolerances(m_rtol, m_atol); + this->setMaxStepSize(m_maxstep); + this->setMaxSteps(m_nmax); + this->setMaxErrTestFails(m_maxErrTestFails); m_integ->initialize(t0, *this); } void ImplicitSurfChem::integrate(doublereal t0, doublereal t1) { - m_integ->initialize(t0, *this); - m_integ->setMaxStepSize(t1 - t0); + this->initialize(t0); + if (fabs(t1 - t0) < m_maxstep || m_maxstep == 0) { + // limit max step size on this run to t1 - t0 + m_integ->setMaxStepSize(t1 - t0); + } m_integ->integrate(t1); updateState(m_integ->solution()); } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 117661e76..d8fc2ea73 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -773,10 +773,12 @@ void InterfaceKinetics::advanceCoverages(doublereal tstep, doublereal rtol, { if (m_integrator == 0) { vector k{this}; - m_integrator = new ImplicitSurfChem(k, rtol, atol, maxStepSize, maxSteps, - maxErrTestFails); - m_integrator->initialize(); + m_integrator = new ImplicitSurfChem(k); } + m_integrator->setTolerances(rtol, atol); + m_integrator->setMaxStepSize(maxStepSize); + m_integrator->setMaxSteps(maxSteps); + m_integrator->setMaxErrTestFails(maxErrTestFails); m_integrator->integrate(0.0, tstep); delete m_integrator; m_integrator = 0;