diff --git a/include/cantera/clib/ctreactor.h b/include/cantera/clib/ctreactor.h
index d535c9913c6..83eac79b8dd 100644
--- a/include/cantera/clib/ctreactor.h
+++ b/include/cantera/clib/ctreactor.h
@@ -53,7 +53,8 @@ extern "C" {
     CANTERA_CAPI int flowdev_del(int i);
     CANTERA_CAPI int flowdev_install(int i, int n, int m);
     CANTERA_CAPI int flowdev_setMaster(int i, int n);
-    CANTERA_CAPI double flowdev_massFlowRate(int i, double time);
+    CANTERA_CAPI double flowdev_massFlowRate2(int i);
+    CANTERA_CAPI double flowdev_massFlowRate(int i, double time); //!< @deprecated To be changed after Cantera 2.5.
     CANTERA_CAPI int flowdev_setMassFlowRate(int i, double mdot);
     CANTERA_CAPI int flowdev_setParameters(int i, int n, const double* v);  //!< @deprecated To be removed after Cantera 2.5.
     CANTERA_CAPI int flowdev_setMassFlowCoeff(int i, double v);
diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h
index 5d84ba4cdb2..46902db7706 100644
--- a/include/cantera/zeroD/FlowDevice.h
+++ b/include/cantera/zeroD/FlowDevice.h
@@ -9,6 +9,7 @@
 #include "cantera/base/ct_defs.h"
 #include "cantera/base/global.h"
 #include "cantera/base/stringUtils.h"
+#include "cantera/base/ctexceptions.h"
 
 namespace Cantera
 {
@@ -54,11 +55,22 @@ class FlowDevice
     }
 
     //! Mass flow rate (kg/s).
+    //! @deprecated The 'time' argument will be removed after Cantera 2.5.
+    //!     Evaluating the mass flow rate at times other than the current time
+    //!     for the reactor network may lead to subtly incorrect results.
     double massFlowRate(double time = -999.0) {
         if (time != -999.0) {
+            warn_deprecated("FlowDevice::massFlowRate", "The 'time' argument"
+                " is deprecated and will be removed after Cantera 2.5.");
             updateMassFlowRate(time);
         }
-        return m_mdot;
+
+        if (m_mdot == Undef) {
+            throw CanteraError("FlowDevice::massFlowRate",
+                               "Flow device is not ready. Try initializing the reactor network.");
+        } else {
+            return m_mdot;
+        }
     }
 
     //! Update the mass flow rate at time 'time'. This must be overloaded in
diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h
index 21720955b98..41528862b28 100644
--- a/include/cantera/zeroD/Reactor.h
+++ b/include/cantera/zeroD/Reactor.h
@@ -184,12 +184,6 @@ class Reactor : public ReactorBase
     //! @param t     the current time
     virtual void evalWalls(double t);
 
-    //! Evaluate inlet and outlet mass flow rates. This is called in evalEqs()
-    //! before setting the state of #m_thermo, since calling the mass flow rate
-    //! functions may modify ThermoPhase objects that are shared with other
-    //! reactors.
-    virtual void evalFlowDevices(double t);
-
     //! Evaluate terms related to surface reactions. Calculates #m_sdot and rate
     //! of change in surface species coverages.
     //! @param t          the current time
@@ -200,6 +194,14 @@ class Reactor : public ReactorBase
     //! Update the state of SurfPhase objects attached to this reactor
     virtual void updateSurfaceState(double* y);
 
+    //! Update the state information needed by connected reactors and flow
+    //! devices. Called from updateState().
+    //! @param updatePressure  Indicates whether to update #m_pressure. Should
+    //!     `true` for reactors where the pressure is a dependent property,
+    //!     calculated from the state, and `false` when the pressure is constant
+    //!     or an independent variable.
+    virtual void updateConnected(bool updatePressure);
+
     //! Get initial conditions for SurfPhase objects attached to this reactor
     virtual void getSurfaceInitialConditions(double* y);
 
diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h
index 2b29a877506..9ac35a5398e 100644
--- a/include/cantera/zeroD/ReactorBase.h
+++ b/include/cantera/zeroD/ReactorBase.h
@@ -269,12 +269,6 @@ class ReactorBase
     vector_fp m_state;
     std::vector<FlowDevice*> m_inlet, m_outlet;
 
-    //! Temporary storage for mass flow rates from each inlet FlowDevice
-    vector_fp m_mdot_in;
-
-    //! Temporary storage for mass flow rates from each outlet FlowDevice
-    vector_fp m_mdot_out;
-
     std::vector<WallBase*> m_wall;
     std::vector<ReactorSurface*> m_surfaces;
     vector_int m_lr;
diff --git a/include/cantera/zeroD/flowControllers.h b/include/cantera/zeroD/flowControllers.h
index a1048de61b6..4d81498bca0 100644
--- a/include/cantera/zeroD/flowControllers.h
+++ b/include/cantera/zeroD/flowControllers.h
@@ -25,10 +25,6 @@ class MassFlowController : public FlowDevice
         return "MassFlowController";
     }
 
-    virtual bool ready() {
-        return FlowDevice::ready() && m_mdot >= 0.0;
-    }
-
     //! Set the fixed mass flow rate (kg/s) through the mass flow controller.
     void setMassFlowRate(double mdot);
 
diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd
index 6dbf992d551..8cc170405a2 100644
--- a/interfaces/cython/cantera/_cantera.pxd
+++ b/interfaces/cython/cantera/_cantera.pxd
@@ -613,6 +613,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
     cdef cppclass CxxFlowDevice "Cantera::FlowDevice":
         CxxFlowDevice()
         string typeStr()
+        double massFlowRate() except +translate_exception
         double massFlowRate(double) except +translate_exception
         cbool install(CxxReactorBase&, CxxReactorBase&) except +translate_exception
         void setPressureFunction(CxxFunc1*) except +translate_exception
diff --git a/interfaces/cython/cantera/examples/reactors/ic_engine.py b/interfaces/cython/cantera/examples/reactors/ic_engine.py
index 8da6a3c5652..13b75c8fefb 100644
--- a/interfaces/cython/cantera/examples/reactors/ic_engine.py
+++ b/interfaces/cython/cantera/examples/reactors/ic_engine.py
@@ -174,8 +174,8 @@ def piston_speed(t):
     states.append(cyl.thermo.state,
                   t=sim.time, ca=crank_angle(sim.time),
                   V=cyl.volume, m=cyl.mass,
-                  mdot_in=inlet_valve.mdot(sim.time),
-                  mdot_out=outlet_valve.mdot(sim.time),
+                  mdot_in=inlet_valve.mass_flow_rate,
+                  mdot_out=outlet_valve.mass_flow_rate,
                   dWv_dt=dWv_dt)
 
 
diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
index 7c16218fd9c..2bde1067d72 100644
--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -707,10 +707,26 @@ cdef class FlowDevice:
         self._upstream = upstream
         self._downstream = downstream
 
-    def mdot(self, double t):
+    property mass_flow_rate:
+        def __get__(self):
+            """
+            The mass flow rate [kg/s] through this device at the current reactor
+            network time.
+            """
+            return self.dev.massFlowRate()
+
+    def mdot(self, double t=-999):
         """
         The mass flow rate [kg/s] through this device at time *t* [s].
+
+        .. deprecated:: 2.5
+
+             To be removed after Cantera 2.5. Replaced with the
+             `mass_flow_rate` property.
         """
+        warnings.warn("To be removed after Cantera 2.5. "
+                "Replaced by property 'mass_flow_rate'", DeprecationWarning)
+
         return self.dev.massFlowRate(t)
 
     def set_pressure_function(self, k):
@@ -799,7 +815,8 @@ cdef class MassFlowController(FlowDevice):
     property mass_flow_rate:
         r"""
         Set the mass flow rate [kg/s] through this controller to be either
-        a constant or an arbitrary function of time. See `Func1`.
+        a constant or an arbitrary function of time. See `Func1`, or get its
+        current value.
 
         Note that depending on the argument type, this method either changes
         the property `mass_flow_coeff` or calls the `set_time_function` method.
@@ -807,6 +824,9 @@ cdef class MassFlowController(FlowDevice):
         >>> mfc.mass_flow_rate = 0.3
         >>> mfc.mass_flow_rate = lambda t: 2.5 * exp(-10 * (t - 0.5)**2)
         """
+        def __get__(self):
+            return self.dev.massFlowRate()
+
         def __set__(self, m):
             if isinstance(m, _numbers.Real):
                 (<CxxMassFlowController*>self.dev).setMassFlowRate(m)
diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py
index 3ce968f0c3a..f7b648bc964 100644
--- a/interfaces/cython/cantera/test/test_reactor.py
+++ b/interfaces/cython/cantera/test/test_reactor.py
@@ -468,13 +468,18 @@ def test_mass_flow_controller(self):
         ma = self.r1.volume * self.r1.density
         Ya = self.r1.Y
 
-        self.assertNear(mfc.mdot(0.1), 0.)
-        self.assertNear(mfc.mdot(0.2), 0.1)
-        self.assertNear(mfc.mdot(1.1), 0.1)
-        self.assertNear(mfc.mdot(1.2), 0.)
-
         self.net.rtol = 1e-11
         self.net.max_time_step = 0.05
+
+        self.net.advance(0.1)
+        self.assertNear(mfc.mass_flow_rate, 0.)
+        self.net.advance(0.2)
+        self.assertNear(mfc.mass_flow_rate, 0.1)
+        self.net.advance(1.1)
+        self.assertNear(mfc.mass_flow_rate, 0.1)
+        self.net.advance(1.2)
+        self.assertNear(mfc.mass_flow_rate, 0.)
+
         self.net.advance(2.5)
 
         mb = self.r1.volume * self.r1.density
@@ -506,8 +511,9 @@ def test_valve1(self):
         self.assertEqual(valve.valve_coeff, k)
         self.assertTrue(self.r1.energy_enabled)
         self.assertTrue(self.r2.energy_enabled)
+        self.net.initialize()
         self.assertNear((self.r1.thermo.P - self.r2.thermo.P) * k,
-                        valve.mdot(0))
+                        valve.mass_flow_rate)
 
         m1a = self.r1.thermo.density * self.r1.volume
         m2a = self.r2.thermo.density * self.r2.volume
@@ -520,7 +526,7 @@ def test_valve1(self):
         m2b = self.r2.thermo.density * self.r2.volume
 
         self.assertNear((self.r1.thermo.P - self.r2.thermo.P) * k,
-                        valve.mdot(0.1))
+                        valve.mass_flow_rate)
         self.assertNear(m1a+m2a, m1b+m2b)
         Y1b = self.r1.thermo.Y
         Y2b = self.r2.thermo.Y
@@ -586,7 +592,7 @@ def speciesMass(k):
             t = self.net.step()
             p1 = self.r1.thermo.P
             p2 = self.r2.thermo.P
-            self.assertNear(mdot(p1-p2), valve.mdot(t))
+            self.assertNear(mdot(p1-p2), valve.mass_flow_rate)
             self.assertArrayNear(Y1, self.r1.Y)
             self.assertNear(speciesMass(kAr), mAr)
             self.assertNear(speciesMass(kO2), mO2)
@@ -600,11 +606,15 @@ def test_valve_timing(self):
         valve.valve_coeff = k
         valve.set_time_function(lambda t: t>.01)
 
-        mdot = valve.valve_coeff * (self.r1.thermo.P - self.r2.thermo.P)
-        self.assertTrue(valve.mdot(0.0)==0.)
-        self.assertTrue(valve.mdot(0.01)==0.)
-        self.assertTrue(valve.mdot(0.01 + 1e-9)==mdot)
-        self.assertTrue(valve.mdot(0.02)==mdot)
+        mdot = lambda: valve.valve_coeff * (self.r1.thermo.P - self.r2.thermo.P)
+        self.net.initialize()
+        self.assertEqual(valve.mass_flow_rate, 0.0)
+        self.net.advance(0.01)
+        self.assertEqual(valve.mass_flow_rate, 0.0)
+        self.net.advance(0.01 + 1e-9)
+        self.assertNear(valve.mass_flow_rate, mdot())
+        self.net.advance(0.02)
+        self.assertNear(valve.mass_flow_rate, mdot())
 
     def test_valve_deprecations(self):
         # Make sure Python deprecation warnings actually get displayed
@@ -669,9 +679,9 @@ def test_pressure_controller1(self):
         t = 0
         while t < 1.0:
             t = self.net.step()
-            self.assertNear(mdot(t), mfc.mdot(t))
+            self.assertNear(mdot(t), mfc.mass_flow_rate)
             dP = self.r1.thermo.P - outlet_reservoir.thermo.P
-            self.assertNear(mdot(t) + 1e-5 * dP, pc.mdot(t))
+            self.assertNear(mdot(t) + 1e-5 * dP, pc.mass_flow_rate)
 
     def test_pressure_controller2(self):
         self.make_reactors(n_reactors=1)
@@ -695,9 +705,9 @@ def test_pressure_controller2(self):
         t = 0
         while t < 1.0:
             t = self.net.step()
-            self.assertNear(mdot(t), mfc.mdot(t))
+            self.assertNear(mdot(t), mfc.mass_flow_rate)
             dP = self.r1.thermo.P - outlet_reservoir.thermo.P
-            self.assertNear(mdot(t) + pfunc(dP), pc.mdot(t))
+            self.assertNear(mdot(t) + pfunc(dP), pc.mass_flow_rate)
 
     def test_pressure_controller_deprecations(self):
         # Make sure Python deprecation warnings actually get displayed
@@ -726,13 +736,13 @@ def test_pressure_controller_errors(self):
 
         p = ct.PressureController(self.r1, self.r2, master=mfc, K=0.5)
 
-        with self.assertRaisesRegex(ct.CanteraError, 'Device is not ready'):
+        with self.assertRaisesRegex(ct.CanteraError, 'is not ready'):
             p = ct.PressureController(self.r1, self.r2, K=0.5)
-            p.mdot(0.0)
+            p.mass_flow_rate
 
-        with self.assertRaisesRegex(ct.CanteraError, 'Device is not ready'):
+        with self.assertRaisesRegex(ct.CanteraError, 'is not ready'):
             p = ct.PressureController(self.r1, self.r2)
-            p.mdot(0.0)
+            p.mass_flow_rate
 
         with self.assertRaisesRegex(ct.CanteraError, 'NotImplementedError'):
             p = ct.PressureController(self.r1, self.r2)
diff --git a/interfaces/matlab/toolbox/@FlowDevice/massFlowRate.m b/interfaces/matlab/toolbox/@FlowDevice/massFlowRate.m
index 8a6d1f6be72..dacd69ca899 100644
--- a/interfaces/matlab/toolbox/@FlowDevice/massFlowRate.m
+++ b/interfaces/matlab/toolbox/@FlowDevice/massFlowRate.m
@@ -9,4 +9,10 @@
 %     The mass flow rate through the :mat:func:`FlowDevice` at the given time
 %
 
-mdot = flowdevicemethods(21, f.index, time);
+if nargin == 1
+    mdot = flowdevicemethods(21, f.index);
+else
+    warning(['"time" argument to massFlowRate is deprecated and will be' ...
+            ' removed after Cantera 2.5.'])
+    mdot = flowdevicemethods(21, f.index, time);
+end
diff --git a/src/clib/ctreactor.cpp b/src/clib/ctreactor.cpp
index 7333f3b2647..13df44541b5 100644
--- a/src/clib/ctreactor.cpp
+++ b/src/clib/ctreactor.cpp
@@ -413,12 +413,24 @@ extern "C" {
     double flowdev_massFlowRate(int i, double time)
     {
         try {
+            warn_deprecated("flowdev_massFlowRate(int i, double time)",
+                "To be changed after Cantera 2.5. 'time' argument will be "
+                "removed. Use flowdev_massFlowRate2(int i) during transition.");
             return FlowDeviceCabinet::item(i).massFlowRate(time);
         } catch (...) {
             return handleAllExceptions(DERR, DERR);
         }
     }
 
+    double flowdev_massFlowRate2(int i)
+    {
+        try {
+            return FlowDeviceCabinet::item(i).massFlowRate();
+        } catch (...) {
+            return handleAllExceptions(DERR, DERR);
+        }
+    }
+
     int flowdev_setMassFlowRate(int i, double mdot)
     {
         /* @deprecated To be removed after Cantera 2.5. */
diff --git a/src/matlab/flowdevicemethods.cpp b/src/matlab/flowdevicemethods.cpp
index fbcb9a43006..d455870321d 100644
--- a/src/matlab/flowdevicemethods.cpp
+++ b/src/matlab/flowdevicemethods.cpp
@@ -79,7 +79,11 @@ void flowdevicemethods(int nlhs, mxArray* plhs[],
         // options that return a value of type 'double'
         switch (job) {
         case 21:
-            r = flowdev_massFlowRate(i, v);
+            if (v == Undef) {
+                r = flowdev_massFlowRate2(i);
+            } else {
+                r = flowdev_massFlowRate(i, v);
+            }
             break;
         default:
             mexErrMsgTxt("unknown job parameter");
diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp
index 5755255d6e3..b2ea8806ced 100644
--- a/src/zeroD/ConstPressureReactor.cpp
+++ b/src/zeroD/ConstPressureReactor.cpp
@@ -55,11 +55,7 @@ void ConstPressureReactor::updateState(doublereal* y)
     }
     m_vol = m_mass / m_thermo->density();
     updateSurfaceState(y + m_nsp + 2);
-
-    // save parameters needed by other connected reactors
-    m_enthalpy = m_thermo->enthalpy_mass();
-    m_intEnergy = m_thermo->intEnergy_mass();
-    m_thermo->saveState(m_state);
+    updateConnected(false);
 }
 
 void ConstPressureReactor::evalEqs(doublereal time, doublereal* y,
@@ -68,7 +64,6 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y,
     double dmdt = 0.0; // dm/dt (gas phase)
     double* dYdt = ydot + 2;
 
-    evalFlowDevices(time);
     evalWalls(time);
     applySensitivity(params);
 
@@ -94,20 +89,22 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y,
     double dHdt = - m_Q;
 
     // add terms for outlets
-    for (size_t i = 0; i < m_outlet.size(); i++) {
-        dmdt -= m_mdot_out[i];
-        dHdt -= m_mdot_out[i] * m_enthalpy;
+    for (auto outlet : m_outlet) {
+        double mdot = outlet->massFlowRate();
+        dmdt -= mdot;
+        dHdt -= mdot * m_enthalpy;
     }
 
     // add terms for inlets
-    for (size_t i = 0; i < m_inlet.size(); i++) {
-        dmdt += m_mdot_in[i]; // mass flow into system
+    for (auto inlet : m_inlet) {
+        double mdot = inlet->massFlowRate();
+        dmdt += mdot; // mass flow into system
         for (size_t n = 0; n < m_nsp; n++) {
-            double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+            double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
             // flow of species into system and dilution by other species
-            dYdt[n] += (mdot_spec - m_mdot_in[i] * Y[n]) / m_mass;
+            dYdt[n] += (mdot_spec - mdot * Y[n]) / m_mass;
         }
-        dHdt += m_mdot_in[i] * m_inlet[i]->enthalpy_mass();
+        dHdt += mdot * inlet->enthalpy_mass();
     }
 
     ydot[0] = dmdt;
diff --git a/src/zeroD/FlowDevice.cpp b/src/zeroD/FlowDevice.cpp
index bd59269e7e8..11ed774d19e 100644
--- a/src/zeroD/FlowDevice.cpp
+++ b/src/zeroD/FlowDevice.cpp
@@ -10,7 +10,7 @@
 namespace Cantera
 {
 
-FlowDevice::FlowDevice() : m_mdot(0.0), m_pfunc(0), m_tfunc(0),
+FlowDevice::FlowDevice() : m_mdot(Undef), m_pfunc(0), m_tfunc(0),
                            m_coeff(1.0), m_type(0),
                            m_nspin(0), m_nspout(0),
                            m_in(0), m_out(0) {}
diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp
index ba6a8ee6869..05db78ac40b 100644
--- a/src/zeroD/IdealGasConstPressureReactor.cpp
+++ b/src/zeroD/IdealGasConstPressureReactor.cpp
@@ -5,6 +5,7 @@
 
 #include "cantera/zeroD/IdealGasConstPressureReactor.h"
 #include "cantera/zeroD/FlowDevice.h"
+#include "cantera/zeroD/ReactorNet.h"
 
 using namespace std;
 
@@ -60,11 +61,7 @@ void IdealGasConstPressureReactor::updateState(doublereal* y)
     m_thermo->setState_TP(y[1], m_pressure);
     m_vol = m_mass / m_thermo->density();
     updateSurfaceState(y + m_nsp + 2);
-
-    // save parameters needed by other connected reactors
-    m_enthalpy = m_thermo->enthalpy_mass();
-    m_intEnergy = m_thermo->intEnergy_mass();
-    m_thermo->saveState(m_state);
+    updateConnected(false);
 }
 
 void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
@@ -74,7 +71,6 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
     double mcpdTdt = 0.0; // m * c_p * dT/dt
     double* dYdt = ydot + 2;
 
-    evalFlowDevices(time);
     applySensitivity(params);
     evalWalls(time);
 
@@ -104,18 +100,19 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
     }
 
     // add terms for outlets
-    for (size_t i = 0; i < m_outlet.size(); i++) {
-        dmdt -= m_mdot_out[i]; // mass flow out of system
+    for (auto outlet : m_outlet) {
+        dmdt -= outlet->massFlowRate(); // mass flow out of system
     }
 
     // add terms for inlets
-    for (size_t i = 0; i < m_inlet.size(); i++) {
-        dmdt += m_mdot_in[i]; // mass flow into system
-        mcpdTdt += m_inlet[i]->enthalpy_mass() * m_mdot_in[i];
+    for (auto inlet : m_inlet) {
+        double mdot = inlet->massFlowRate();
+        dmdt += mdot; // mass flow into system
+        mcpdTdt += inlet->enthalpy_mass() * mdot;
         for (size_t n = 0; n < m_nsp; n++) {
-            double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+            double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
             // flow of species into system and dilution by other species
-            dYdt[n] += (mdot_spec - m_mdot_in[i] * Y[n]) / m_mass;
+            dYdt[n] += (mdot_spec - mdot * Y[n]) / m_mass;
             mcpdTdt -= m_hk[n] / mw[n] * mdot_spec;
         }
     }
diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp
index 20e1af4cad7..93553b0bc29 100644
--- a/src/zeroD/IdealGasReactor.cpp
+++ b/src/zeroD/IdealGasReactor.cpp
@@ -65,12 +65,7 @@ void IdealGasReactor::updateState(doublereal* y)
     m_thermo->setMassFractions_NoNorm(y+3);
     m_thermo->setState_TR(y[2], m_mass / m_vol);
     updateSurfaceState(y + m_nsp + 3);
-
-    // save parameters needed by other connected reactors
-    m_enthalpy = m_thermo->enthalpy_mass();
-    m_pressure = m_thermo->pressure();
-    m_intEnergy = m_thermo->intEnergy_mass();
-    m_thermo->saveState(m_state);
+    updateConnected(true);
 }
 
 void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
@@ -80,7 +75,6 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
     double mcvdTdt = 0.0; // m * c_v * dT/dt
     double* dYdt = ydot + 3;
 
-    evalFlowDevices(time);
     evalWalls(time);
     applySensitivity(params);
     m_thermo->restoreState(m_state);
@@ -109,21 +103,21 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
     }
 
     // add terms for outlets
-    for (size_t i = 0; i < m_outlet.size(); i++) {
-        // double mdot_out = m_outlet[i]->massFlowRate(time);
-        dmdt -= m_mdot_out[i]; // mass flow out of system
-        mcvdTdt -= m_mdot_out[i] * m_pressure * m_vol / m_mass; // flow work
+    for (auto outlet : m_outlet) {
+        double mdot = outlet->massFlowRate();
+        dmdt -= mdot; // mass flow out of system
+        mcvdTdt -= mdot * m_pressure * m_vol / m_mass; // flow work
     }
 
     // add terms for inlets
-    for (size_t i = 0; i < m_inlet.size(); i++) {
-        // double mdot_in = m_inlet[i]->massFlowRate(time);
-        dmdt += m_mdot_in[i]; // mass flow into system
-        mcvdTdt += m_inlet[i]->enthalpy_mass() * m_mdot_in[i];
+    for (auto inlet : m_inlet) {
+        double mdot = inlet->massFlowRate();
+        dmdt += mdot; // mass flow into system
+        mcvdTdt += inlet->enthalpy_mass() * mdot;
         for (size_t n = 0; n < m_nsp; n++) {
-            double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+            double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
             // flow of species into system and dilution by other species
-            dYdt[n] += (mdot_spec - m_mdot_in[i] * Y[n]) / m_mass;
+            dYdt[n] += (mdot_spec - mdot * Y[n]) / m_mass;
 
             // In combination with h_in*mdot_in, flow work plus thermal
             // energy carried with the species
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index 02856b98235..85fd70a3520 100644
--- a/src/zeroD/Reactor.cpp
+++ b/src/zeroD/Reactor.cpp
@@ -81,10 +81,7 @@ void Reactor::initialize(doublereal t0)
     m_thermo->restoreState(m_state);
     m_sdot.resize(m_nsp, 0.0);
     m_wdot.resize(m_nsp, 0.0);
-
-    m_enthalpy = m_thermo->enthalpy_mass();
-    m_pressure = m_thermo->pressure();
-    m_intEnergy = m_thermo->intEnergy_mass();
+    updateConnected(true);
 
     for (size_t n = 0; n < m_wall.size(); n++) {
         WallBase* W = m_wall[n];
@@ -165,12 +162,7 @@ void Reactor::updateState(doublereal* y)
     }
 
     updateSurfaceState(y + m_nsp + 3);
-
-    // save parameters needed by other connected reactors
-    m_enthalpy = m_thermo->enthalpy_mass();
-    m_pressure = m_thermo->pressure();
-    m_intEnergy = m_thermo->intEnergy_mass();
-    m_thermo->saveState(m_state);
+    updateConnected(true);
 }
 
 void Reactor::updateSurfaceState(double* y)
@@ -182,13 +174,31 @@ void Reactor::updateSurfaceState(double* y)
     }
 }
 
+void Reactor::updateConnected(bool updatePressure) {
+    // save parameters needed by other connected reactors
+    m_enthalpy = m_thermo->enthalpy_mass();
+    if (updatePressure) {
+        m_pressure = m_thermo->pressure();
+    }
+    m_intEnergy = m_thermo->intEnergy_mass();
+    m_thermo->saveState(m_state);
+
+    // Update the mass flow rate of connected flow devices
+    double time = m_net->time();
+    for (size_t i = 0; i < m_outlet.size(); i++) {
+        m_outlet[i]->updateMassFlowRate(time);
+    }
+    for (size_t i = 0; i < m_inlet.size(); i++) {
+        m_inlet[i]->updateMassFlowRate(time);
+    }
+}
+
 void Reactor::evalEqs(doublereal time, doublereal* y,
                       doublereal* ydot, doublereal* params)
 {
     double dmdt = 0.0; // dm/dt (gas phase)
     double* dYdt = ydot + 3;
 
-    evalFlowDevices(time);
     evalWalls(time);
     applySensitivity(params);
     m_thermo->restoreState(m_state);
@@ -223,23 +233,25 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
     }
 
     // add terms for outlets
-    for (size_t i = 0; i < m_outlet.size(); i++) {
-        dmdt -= m_mdot_out[i]; // mass flow out of system
+    for (auto outlet : m_outlet) {
+        double mdot = outlet->massFlowRate();
+        dmdt -= mdot; // mass flow out of system
         if (m_energy) {
-            ydot[2] -= m_mdot_out[i] * m_enthalpy;
+            ydot[2] -= mdot * m_enthalpy;
         }
     }
 
     // add terms for inlets
-    for (size_t i = 0; i < m_inlet.size(); i++) {
-        dmdt += m_mdot_in[i]; // mass flow into system
+    for (auto inlet : m_inlet) {
+        double mdot = inlet->massFlowRate();
+        dmdt += mdot; // mass flow into system
         for (size_t n = 0; n < m_nsp; n++) {
-            double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
+            double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
             // flow of species into system and dilution by other species
-            dYdt[n] += (mdot_spec - m_mdot_in[i] * Y[n]) / m_mass;
+            dYdt[n] += (mdot_spec - mdot * Y[n]) / m_mass;
         }
         if (m_energy) {
-            ydot[2] += m_mdot_in[i] * m_inlet[i]->enthalpy_mass();
+            ydot[2] += mdot * inlet->enthalpy_mass();
         }
     }
 
@@ -258,16 +270,6 @@ void Reactor::evalWalls(double t)
     }
 }
 
-void Reactor::evalFlowDevices(double t)
-{
-    for (size_t i = 0; i < m_outlet.size(); i++) {
-        m_mdot_out[i] = m_outlet[i]->massFlowRate(t);
-    }
-    for (size_t i = 0; i < m_inlet.size(); i++) {
-        m_mdot_in[i] = m_inlet[i]->massFlowRate(t);
-    }
-}
-
 double Reactor::evalSurfaces(double t, double* ydot)
 {
     const vector_fp& mw = m_thermo->molecularWeights();
diff --git a/src/zeroD/ReactorBase.cpp b/src/zeroD/ReactorBase.cpp
index 566250e457d..edcb7e1532c 100644
--- a/src/zeroD/ReactorBase.cpp
+++ b/src/zeroD/ReactorBase.cpp
@@ -47,13 +47,11 @@ void ReactorBase::syncState()
 void ReactorBase::addInlet(FlowDevice& inlet)
 {
     m_inlet.push_back(&inlet);
-    m_mdot_in.push_back(0.0);
 }
 
 void ReactorBase::addOutlet(FlowDevice& outlet)
 {
     m_outlet.push_back(&outlet);
-    m_mdot_out.push_back(0.0);
 }
 
 void ReactorBase::addWall(WallBase& w, int lr)
diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp
index d90d819ccf4..e51627c7121 100644
--- a/src/zeroD/ReactorNet.cpp
+++ b/src/zeroD/ReactorNet.cpp
@@ -233,6 +233,7 @@ void ReactorNet::addReactor(Reactor& r)
 void ReactorNet::eval(doublereal t, doublereal* y,
                       doublereal* ydot, doublereal* p)
 {
+    m_time = t; // This will be replaced at the end of the timestep
     updateState(y);
     for (size_t n = 0; n < m_reactors.size(); n++) {
         m_reactors[n]->evalEqs(t, y + m_start[n], ydot + m_start[n], p);
diff --git a/src/zeroD/flowControllers.cpp b/src/zeroD/flowControllers.cpp
index 35095658cb3..f4018819510 100644
--- a/src/zeroD/flowControllers.cpp
+++ b/src/zeroD/flowControllers.cpp
@@ -52,7 +52,8 @@ void PressureController::updateMassFlowRate(double time)
     } else {
         mdot *= delta_P;
     }
-    mdot += m_master->massFlowRate(time);
+    m_master->updateMassFlowRate(time);
+    mdot += m_master->massFlowRate();
     m_mdot = std::max(mdot, 0.0);
 }