[Reactor] Instantiate Integrator early to simplify setting options

include/cantera/zeroD/FlowReactor.h

@@ -22,6 +22,10 @@ class FlowReactor : public IdealGasReactor
         return "FlowReactor";
     }
 
+    bool isOde() const override {
+        return false;
+    }
+
     void getState(double* y) override {
         throw NotImplementedError("FlowReactor::getState");
     }

include/cantera/zeroD/Reactor.h

@@ -49,6 +49,13 @@ class Reactor : public ReactorBase
         return "Reactor";
     }
 
+    //! Indicate whether the governing equations for this reactor type are a system of
+    //! ODEs or DAEs. In the first case, this class implements the eval() method. In the
+    //! second case, this class implements the evalDae() method.
+    virtual bool isOde() const {
+        return true;
+    }
+
     /**
      * Insert something into the reactor. The 'something' must belong to a class
      * that is a subclass of both ThermoPhase and Kinetics.

include/cantera/zeroD/ReactorNet.h

@@ -143,10 +143,9 @@ class ReactorNet : public FuncEval
         suppressErrors(!m_verbose);
     }
 
-    //! Return a reference to the integrator.
-    Integrator& integrator() {
-        return *m_integ;
-    }
+    //! Return a reference to the integrator. Only valid after adding at least one
+    //! reactor to the network.
+    Integrator& integrator();
 
     //! Update the state of all the reactors in the network to correspond to
     //! the values in the solution vector *y*.
@@ -320,17 +319,12 @@ class ReactorNet : public FuncEval
     double m_atolsens = 1.0e-6;
     shared_ptr<PreconditionerBase> m_precon;
     string m_linearSolverType;
-    int m_maxSteps = -1;
 
     //! Maximum integrator internal timestep. Default of 0.0 means infinity.
     double m_maxstep = 0.0;
 
-    int m_maxErrTestFails = 0;
     bool m_verbose = false;
 
-    // flag indicating whether this is a flow reactor net or a regular reactor net
-    bool m_is_dae;
-
     //! Names corresponding to each sensitivity parameter
     std::vector<std::string> m_paramNames;
 

interfaces/cython/cantera/reactor.pxd

@@ -137,7 +137,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
 
     cdef cppclass CxxReactorNet "Cantera::ReactorNet":
         CxxReactorNet()
-        void addReactor(CxxReactor&)
+        void addReactor(CxxReactor&) except +translate_exception
         double advance(double, cbool) except +translate_exception
         double step() except +translate_exception
         void initialize() except +translate_exception
@@ -148,9 +148,9 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         double rtol()
         double atol()
         double maxTimeStep()
-        void setMaxTimeStep(double)
-        void setMaxErrTestFails(int)
-        void setMaxSteps(int)
+        void setMaxTimeStep(double) except +translate_exception
+        void setMaxErrTestFails(int) except +translate_exception
+        void setMaxSteps(int) except +translate_exception
         int maxSteps()
         cbool verbose()
         void setVerbose(cbool)

src/zeroD/ReactorNet.cpp

@@ -16,8 +16,7 @@
 namespace Cantera
 {
 
-ReactorNet::ReactorNet() :
-    m_is_dae(false)
+ReactorNet::ReactorNet()
 {
     suppressErrors(true);
 }
@@ -35,13 +34,12 @@ void ReactorNet::setInitialTime(double time)
 void ReactorNet::setMaxTimeStep(double maxstep)
 {
     m_maxstep = maxstep;
-    m_init = false;
+    integrator().setMaxStepSize(m_maxstep);
 }
 
 void ReactorNet::setMaxErrTestFails(int nmax)
 {
-    m_maxErrTestFails = nmax;
-    m_init = false;
+    integrator().setMaxErrTestFails(nmax);
 }
 
 void ReactorNet::setTolerances(double rtol, double atol)
@@ -69,12 +67,6 @@ void ReactorNet::setSensitivityTolerances(double rtol, double atol)
 void ReactorNet::initialize()
 {
     m_nv = 0;
-    m_integ.reset(newIntegrator(m_is_dae ? "IDA" : "CVODE"));
-    // use backward differencing, with a full Jacobian computed
-    // numerically, and use a Newton linear iterator
-    m_integ->setMethod(BDF_Method);
-    m_integ->setLinearSolverType("DENSE");
-
     debuglog("Initializing reactor network.\n", m_verbose);
     if (m_reactors.empty()) {
         throw CanteraError("ReactorNet::initialize",
@@ -105,8 +97,6 @@ void ReactorNet::initialize()
     fill(m_atol.begin(), m_atol.end(), m_atols);
     m_integ->setTolerances(m_rtol, neq(), m_atol.data());
     m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens);
-    m_integ->setMaxStepSize(m_maxstep);
-    m_integ->setMaxErrTestFails(m_maxErrTestFails);
     if (!m_linearSolverType.empty()) {
         m_integ->setLinearSolverType(m_linearSolverType);
     }
@@ -121,9 +111,6 @@ void ReactorNet::initialize()
     if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) {
         checkPreconditionerSupported();
     }
-    if (m_maxSteps != -1) {
-        m_integ->setMaxSteps(m_maxSteps);
-    }
     m_integrator_init = true;
     m_init = true;
 }
@@ -156,19 +143,12 @@ void ReactorNet::setPreconditioner(shared_ptr<PreconditionerBase> preconditioner
 
 void ReactorNet::setMaxSteps(int nmax)
 {
-    m_maxSteps = nmax;
-    if (m_integ) {
-        m_integ->setMaxSteps(nmax);
-    }
+    integrator().setMaxSteps(nmax);
 }
 
 int ReactorNet::maxSteps()
 {
-    if (m_integ) {
-        return m_integ->maxSteps();
-    } else {
-        return m_maxSteps;
-    }
+    return integrator().maxSteps();
 }
 
 void ReactorNet::advance(doublereal time)
@@ -284,19 +264,30 @@ int ReactorNet::lastOrder()
 
 void ReactorNet::addReactor(Reactor& r)
 {
-    r.setNetwork(this);
-    try {
-        dynamic_cast<FlowReactor&>(r);
-        m_is_dae = true;
-    } 
-    catch(std::bad_cast&) {
-        if (m_is_dae)
-        {
+    for (auto current : m_reactors) {
+        if (current->isOde() != r.isOde()) {
             throw CanteraError("ReactorNet::addReactor",
-                               "Cannot mix Reactors and FlowReactor's");
+                "Cannot mix Reactor types using both ODEs and DAEs ({} and {})",
+                current->type(), r.type());
         }
     }
+    r.setNetwork(this);
     m_reactors.push_back(&r);
+    if (!m_integ) {
+        m_integ.reset(newIntegrator(r.isOde() ? "CVODE" : "IDA"));
+        // use backward differencing, with a full Jacobian computed
+        // numerically, and use a Newton linear iterator
+        m_integ->setMethod(BDF_Method);
+        m_integ->setLinearSolverType("DENSE");
+    }
+}
+
+Integrator& ReactorNet::integrator() {
+    if (m_integ == nullptr) {
+        throw CanteraError("ReactorNet::integrator",
+            "Integrator has not been instantiated. Add one or more reactors first.");
+    }
+    return *m_integ;
 }
 
 void ReactorNet::eval(double t, double* y, double* ydot, double* p)

test/python/test_reactor.py

@@ -1513,6 +1513,13 @@ def make_reactors(self, gas, surf):
         sim = ct.ReactorNet([r])
         return r, rsurf, sim
 
+    def test_mixed_reactor_types(self):
+        surf, gas = self.import_phases()
+        r1 = ct.FlowReactor(gas)
+        r2 = ct.IdealGasReactor(gas)
+        with pytest.raises(ct.CanteraError, match="Cannot mix Reactor types"):
+            ct.ReactorNet([r1, r2])
+
     def test_unrecoverable_integrator_errors(self):
         surf, gas = self.import_phases()