[Reactor] Distinguish between time and distance as independent variables

include/cantera/zeroD/FlowReactor.h

@@ -25,6 +25,10 @@ class FlowReactor : public IdealGasReactor
         return false;
     }
 
+    bool timeIsIndependent() const override {
+        return false;
+    }
+
     //! Not implemented; FlowReactor implements getStateDAE() instead.
     void getState(double* y) override {
         throw NotImplementedError("FlowReactor::getState");

include/cantera/zeroD/Reactor.h

@@ -56,6 +56,12 @@ class Reactor : public ReactorBase
         return true;
     }
 
+    //! Indicates whether the governing equations for this reactor are functions of time
+    //! or a spatial variable. All reactors in a network must have the same value.
+    virtual bool timeIsIndependent() 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

@@ -19,9 +19,10 @@ class PreconditionerBase;
 
 //! A class representing a network of connected reactors.
 /*!
- *  This class is used to integrate the time-dependent governing equations for
- *  a network of reactors (Reactor, ConstPressureReactor) connected by various
- *  means, for example Wall, MassFlowController, Valve, or PressureController.
+ *  This class is used to integrate the governing equations for a network of reactors
+ *  that are time dependent (Reactor, ConstPressureReactor) connected by various
+ *  means, for example Wall, MassFlowController, Valve, or PressureController; or
+ *  reactors dependent on a single spatial variable (FlowReactor).
  *
  * @ingroup ZeroD
  */
@@ -45,20 +46,21 @@ class ReactorNet : public FuncEval
     //! @param preconditioner preconditioner object used for the linear solver
     void setPreconditioner(shared_ptr<PreconditionerBase> preconditioner);
 
-    //! Set initial time. Default = 0.0 s. Restarts integration from this time
-    //! using the current mixture state as the initial condition.
+    //! Set the initial value of the independent variable (typically time).
+    //! Default = 0.0 s. Restarts integration from this value using the current mixture
+    //! state as the initial condition.
     void setInitialTime(double time);
 
-    //! Get the maximum time step.
+    //! Get the maximum integrator step.
     double maxTimeStep() {
         return m_maxstep;
     }
 
-    //! Set the maximum time step.
+    //! Set the maximum integrator step.
     void setMaxTimeStep(double maxstep);
 
     //! Set the maximum number of error test failures permitted by the CVODES
-    //! integrator in a single time step.
+    //! integrator in a single step.
     void setMaxErrTestFails(int nmax);
 
     //! Set the relative and absolute tolerances for the integrator.
@@ -68,10 +70,13 @@ class ReactorNet : public FuncEval
     //! sensitivity equations.
     void setSensitivityTolerances(double rtol, double atol);
 
-    //! Current value of the simulation time.
-    double time() {
-        return m_time;
-    }
+    //! Current value of the simulation time [s], for reactor networks that are solved
+    //! in the time domain.
+    double time();
+
+    //! Current position [m] along the length of the reactor network, for reactors that
+    //! are solved as a function of space.
+    double distance();
 
     //! Relative tolerance.
     double rtol() {
@@ -96,29 +101,30 @@ class ReactorNet : public FuncEval
     //! Problem type of integrator
     std::string linearSolverType() const;
 
-    //! Returns the maximum number of internal integration time-steps the
-    //!  integrator will take before reaching the next output time
+    //! Returns the maximum number of internal integration steps the
+    //! integrator will take before reaching the next output point
     int maxSteps();
 
     /**
-     * Advance the state of all reactors in time. Take as many internal
-     * timesteps as necessary to reach *time*.
-     * @param time Time to advance to (s).
+     * Advance the state of all reactors in the independent variable (time or space).
+     * Take as many internal steps as necessary to reach *t*.
+     * @param t Time/distance to advance to (s or m).
      */
-    void advance(double time);
+    void advance(double t);
 
     /**
-     * Advance the state of all reactors in time. Take as many internal
-     * timesteps as necessary towards *time*. If *applylimit* is true,
-     * the advance step will be automatically reduced if needed to
-     * stay within limits (set by setAdvanceLimit).
-     * Returns the time at the end of integration.
-     * @param time Time to advance to (s).
+     * Advance the state of all reactors in the independent variable (time or space).
+     * Take as many internal steps as necessary towards *t*. If *applylimit* is true,
+     * the advance step will be automatically reduced if needed to stay within limits
+     * (set by setAdvanceLimit).
+     * Returns the time/distance at the end of integration.
+     * @param t Time/distance to advance to (s or m).
      * @param applylimit Limit advance step (boolean).
      */
-    double advance(double time, bool applylimit);
+    double advance(double t, bool applylimit);
 
-    //! Advance the state of all reactors in time.
+    //! Advance the state of all reactors with respect to the independent variable
+    //! (time or space). Returns the new value of the independent variable [s or m].
     double step();
 
     //! Add the reactor *r* to this reactor network.
@@ -180,9 +186,10 @@ class ReactorNet : public FuncEval
 
     //! Evaluate the Jacobian matrix for the reactor network.
     /*!
-     *  @param[in] t Time at which to evaluate the Jacobian
-     *  @param[in] y Global state vector at time *t*
-     *  @param[out] ydot Time derivative of the state vector evaluated at *t*.
+     *  @param[in] t Time/distance at which to evaluate the Jacobian
+     *  @param[in] y Global state vector at *t*
+     *  @param[out] ydot Derivative of the state vector evaluated at *t*, with respect
+     *      to *t*.
      *  @param[in] p sensitivity parameter vector (unused?)
      *  @param[out] j Jacobian matrix, size neq() by neq().
      */
@@ -207,7 +214,7 @@ class ReactorNet : public FuncEval
     virtual void getState(double* y);
     virtual void getStateDae(double* y, double* ydot);
 
-    //! Return k-th derivative at the current time
+    //! Return k-th derivative at the current state of the system
     virtual void getDerivative(int k, double* dky);
 
     virtual void getConstraints(double* constraints);
@@ -259,8 +266,8 @@ class ReactorNet : public FuncEval
         m_integrator_init = false;
     }
 
-    //! Set the maximum number of internal integration time-steps the
-    //! integrator will take before reaching the next output time
+    //! Set the maximum number of internal integration steps the
+    //! integrator will take before reaching the next output point
     //! @param nmax The maximum number of steps, setting this value
     //!             to zero disables this option.
     virtual void setMaxSteps(int nmax);
@@ -304,7 +311,11 @@ class ReactorNet : public FuncEval
 
     std::vector<Reactor*> m_reactors;
     std::unique_ptr<Integrator> m_integ;
+
+    //! The independent variable in the system. May be either time or space depending
+    //! on the type of reactors in the network.
     double m_time = 0.0;
+
     bool m_init = false;
     bool m_integrator_init = false; //!< True if integrator initialization is current
     size_t m_nv = 0;
@@ -325,6 +336,9 @@ class ReactorNet : public FuncEval
 
     bool m_verbose = false;
 
+    //! Indicates whether time or space is the independent variable
+    bool m_timeIsIndependent = true;
+
     //! Names corresponding to each sensitivity parameter
     std::vector<std::string> m_paramNames;
 

interfaces/cython/cantera/reactor.pxd

@@ -161,7 +161,8 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         double step() except +translate_exception
         void initialize() except +translate_exception
         void reinitialize() except +translate_exception
-        double time()
+        double time() except +translate_exception
+        double distance() except +translate_exception
         void setInitialTime(double)
         void setTolerances(double, double)
         double rtol()

interfaces/cython/cantera/reactor.pyx

@@ -440,7 +440,7 @@ cdef class IdealGasConstPressureMoleReactor(Reactor):
 cdef class FlowReactor(Reactor):
     """
     A steady-state plug flow reactor with constant cross sectional area.
-    Time integration follows a fluid element along the length of the reactor.
+    Integration follows a fluid element along the length of the reactor.
     The reactor is assumed to be frictionless and adiabatic.
     """
     reactor_type = "FlowReactor"
@@ -1291,19 +1291,22 @@ cdef class ReactorNet:
 
     def advance(self, double t, pybool apply_limit=True):
         """
-        Advance the state of the reactor network in time from the current time
-        towards time ``t`` in seconds, taking as many integrator time steps as necessary.
-        If ``apply_limit`` is true and an advance limit is specified, the reactor
-        state at the end of the timestep is estimated prior to advancing. If
-        the difference exceed limits, the end time is reduced by half until
-        the projected end state remains within specified limits.
-        Returns the time reached at the end of integration.
+        Advance the state of the reactor network from the current time/distance towards
+        the specified value ``t`` of the independent variable, which depends on the type
+        of reactors included in the network.
+
+        The integrator will take as many steps as necessary to reach ``t``. If
+        ``apply_limit`` is true and an advance limit is specified, the reactor state at
+        the end of the step is estimated prior to advancing. If the difference exceed
+        limits, the end value is reduced by half until the projected end state remains
+        within specified limits. Returns the time/distance reached at the end of
+        integration.
         """
         return self.net.advance(t, apply_limit)
 
     def step(self):
         """
-        Take a single internal time step. The time after taking the step is
+        Take a single internal step. The time/distance after taking the step is
         returned.
         """
         return self.net.step()
@@ -1322,10 +1325,20 @@ cdef class ReactorNet:
         """
         self.net.reinitialize()
 
-    property time:
-        """The current time [s]."""
-        def __get__(self):
-            return self.net.time()
+    @property
+    def time(self):
+        """
+        The current time [s], for reactor networks that are solved in the time domain.
+        """
+        return self.net.time()
+
+    @property
+    def distance(self):
+        """
+        The current distance[ m] along the length of the reactor network, for reactors
+        that are solved as a function of space.
+        """
+        return self.net.distance()
 
     def set_initial_time(self, double t):
         """
@@ -1348,8 +1361,8 @@ cdef class ReactorNet:
 
     property max_err_test_fails:
         """
-        The maximum number of error test failures permitted by the CVODES
-        integrator in a single time step. The default is 10.
+        The maximum number of error test failures permitted by the CVODES integrator
+        in a single step. The default is 10.
         """
         def __set__(self, n):
             self.net.setMaxErrTestFails(n)
@@ -1404,8 +1417,8 @@ cdef class ReactorNet:
 
     property max_steps:
         """
-        The maximum number of internal integration time-steps that CVODES
-        is allowed to take before reaching the next output time.
+        The maximum number of internal integration steps that CVODES
+        is allowed to take before reaching the next output point.
         """
         def __set__(self, nsteps):
             self.net.setMaxSteps(nsteps)
@@ -1486,7 +1499,7 @@ cdef class ReactorNet:
         string or an integer. See `component_index` and `sensitivities` to
         determine the integer index for the variables and the definition of the
         resulting sensitivity coefficient. If it is not given, ``r`` defaults to
-        the first reactor. Returns an empty array until the first time step is
+        the first reactor. Returns an empty array until the first integration step is
         taken.
         """
         if isinstance(component, int):
@@ -1508,7 +1521,7 @@ cdef class ReactorNet:
         reversible reactions).
 
         The sensitivities are returned in an array with dimensions *(n_vars,
-        n_sensitivity_params)*, unless no timesteps have been taken, in which
+        n_sensitivity_params)*, unless no integration steps have been taken, in which
         case the shape is *(0, n_sensitivity_params)*. The order of the
         variables (that is, rows) is:
 
@@ -1578,7 +1591,8 @@ cdef class ReactorNet:
 
     def get_derivative(self, k):
         """
-        Get the k-th time derivative of the state vector of the reactor network.
+        Get the k-th derivative of the state vector of the reactor network with respect
+        to the independent integrator variable (time/distance).
         """
         if not self.n_vars:
             raise CanteraError('ReactorNet empty or not initialized.')

samples/python/reactors/surf_pfr.py

@@ -64,8 +64,8 @@ class and the SUNDIALS IDA solver, in contrast to the approximation as a chain o
 print('  {:10f}  {:10f}  {:10f}  {:10f}'.format(
       0, *r.thermo['CH4', 'H2', 'CO'].X))
 
-while sim.time < length:
-    dist = sim.time * 1e3
+while sim.distance < length:
+    dist = sim.distance * 1e3  # convert to mm
     sim.step()
 
     if n % 100 == 0 or (dist > 1 and n % 10 == 0):

samples/python/surface_chemistry/1D_pfr_surfchem.py

@@ -48,14 +48,13 @@
 kN = gas_si_n_interface.kinetics_species_index('N(D)')
 kSi = gas_si_n_interface.kinetics_species_index('Si(D)')
 
-# Integrate the reactor network. Note that "time" is really the coordinate along
-# the length of the plug flow reactor
-while net.time < 0.6:
-    print(net.time, rsurf.coverages)
+# Integrate the reactor network
+while net.distance < 0.6:
+    print(net.distance, rsurf.coverages)
     net.step()
     wdot = rsurf.kinetics.net_production_rates
     soln.append(TDY=reactor.thermo.TDY,
-                x=net.time,
+                x=net.distance,
                 speed=reactor.speed,
                 surf_coverages=rsurf.coverages,
                 N_dep=wdot[kN],

src/zeroD/Reactor.cpp

@@ -191,7 +191,10 @@ void Reactor::updateConnected(bool updatePressure) {
     m_thermo->saveState(m_state);
 
     // Update the mass flow rate of connected flow devices
-    double time = (m_net != nullptr) ? m_net->time() : 0.0;
+    double time = 0.0;
+    if (m_net != nullptr) {
+        time = (timeIsIndependent()) ? m_net->time() : m_net->distance();
+    }
     for (size_t i = 0; i < m_outlet.size(); i++) {
         m_outlet[i]->updateMassFlowRate(time);
     }

src/zeroD/ReactorNet.cpp

@@ -64,6 +64,24 @@ void ReactorNet::setSensitivityTolerances(double rtol, double atol)
     m_init = false;
 }
 
+double ReactorNet::time() {
+    if (m_timeIsIndependent) {
+        return m_time;
+    } else {
+        throw CanteraError("ReactorNet::time", "Time is not the independent variable"
+            " for this reactor network.");
+    }
+}
+
+double ReactorNet::distance() {
+    if (!m_timeIsIndependent) {
+        return m_time;
+    } else {
+        throw CanteraError("ReactorNet::distance", "Distance is not the independent"
+            " variable for this reactor network.");
+    }
+}
+
 void ReactorNet::initialize()
 {
     m_nv = 0;
@@ -270,7 +288,13 @@ void ReactorNet::addReactor(Reactor& r)
                 "Cannot mix Reactor types using both ODEs and DAEs ({} and {})",
                 current->type(), r.type());
         }
+        if (current->timeIsIndependent() != r.timeIsIndependent()) {
+            throw CanteraError("ReactorNet::addReactor",
+                "Cannot mix Reactor types using time and space as independent variables"
+                "\n({} and {})", current->type(), r.type());
+        }
     }
+    m_timeIsIndependent = r.timeIsIndependent();
     r.setNetwork(this);
     m_reactors.push_back(&r);
     if (!m_integ) {