From 0e3317a5106d0a0e18c13d28acd5bd0ffdbd3d20 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 29 Jun 2020 09:04:18 -0400 Subject: [PATCH 1/3] [Reactor] Update mass flow rates after ReactorNet.step --- include/cantera/zeroD/Reactor.h | 8 ++++++ src/zeroD/ConstPressureReactor.cpp | 6 +---- src/zeroD/IdealGasConstPressureReactor.cpp | 7 ++--- src/zeroD/IdealGasReactor.cpp | 7 +---- src/zeroD/Reactor.cpp | 31 +++++++++++++++------- 5 files changed, 33 insertions(+), 26 deletions(-) diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 21720955b98..0e16238460d 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -200,6 +200,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/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 5755255d6e3..d95ac064114 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, diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index ba6a8ee6869..9bd057f8cd1 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, diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 20e1af4cad7..69cfa1a8073 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, diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 02856b98235..2a94f8b497c 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,6 +174,25 @@ 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 (auto device : m_inlet) { + device->updateMassFlowRate(time); + } + for (auto device : m_outlet) { + device->updateMassFlowRate(time); + } +} + void Reactor::evalEqs(doublereal time, doublereal* y, doublereal* ydot, doublereal* params) { From e835772b6959173f9204321d55940d3bec0ac856 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 29 Jun 2020 16:43:06 -0400 Subject: [PATCH 2/3] [Reactor] Deprecate 'time' argument to FlowDevice::massFlowRate Since the mass flow rate can depend on both a function of time and the state of the connected reactors, the only time for which the mass flow rate can be reliably returned is the current time of the reactor network. --- include/cantera/clib/ctreactor.h | 3 +- include/cantera/zeroD/FlowDevice.h | 14 ++++- include/cantera/zeroD/flowControllers.h | 4 -- interfaces/cython/cantera/_cantera.pxd | 1 + .../cantera/examples/reactors/ic_engine.py | 4 +- interfaces/cython/cantera/reactor.pyx | 24 ++++++++- .../cython/cantera/test/test_reactor.py | 52 +++++++++++-------- .../matlab/toolbox/@FlowDevice/massFlowRate.m | 8 ++- src/clib/ctreactor.cpp | 12 +++++ src/matlab/flowdevicemethods.cpp | 6 ++- src/zeroD/FlowDevice.cpp | 2 +- src/zeroD/Reactor.cpp | 6 ++- src/zeroD/flowControllers.cpp | 3 +- 13 files changed, 102 insertions(+), 37 deletions(-) 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/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): (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/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/Reactor.cpp b/src/zeroD/Reactor.cpp index 2a94f8b497c..fb821adaa03 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -272,10 +272,12 @@ 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); + m_outlet[i]->updateMassFlowRate(t); + m_mdot_out[i] = m_outlet[i]->massFlowRate(); } for (size_t i = 0; i < m_inlet.size(); i++) { - m_mdot_in[i] = m_inlet[i]->massFlowRate(t); + m_inlet[i]->updateMassFlowRate(t); + m_mdot_in[i] = m_inlet[i]->massFlowRate(); } } 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); } From d58fdc38da4f4f6b245eea6d71673591abc5fe91 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 29 Jun 2020 17:27:44 -0400 Subject: [PATCH 3/3] [Reactor] Simplify updates of flow device mass flow rates Eliminate unnecessary variables m_mdot_in and m_mdot_out, and remove the need to update the mass flow rates in both updateState and evalEqs by using the current (trial) time when updating them in updateState. --- include/cantera/zeroD/Reactor.h | 6 ---- include/cantera/zeroD/ReactorBase.h | 6 ---- src/zeroD/ConstPressureReactor.cpp | 19 ++++++----- src/zeroD/IdealGasConstPressureReactor.cpp | 16 ++++----- src/zeroD/IdealGasReactor.cpp | 21 ++++++------ src/zeroD/Reactor.cpp | 39 ++++++++-------------- src/zeroD/ReactorBase.cpp | 2 -- src/zeroD/ReactorNet.cpp | 1 + 8 files changed, 43 insertions(+), 67 deletions(-) diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 0e16238460d..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 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 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 m_wall; std::vector m_surfaces; vector_int m_lr; diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index d95ac064114..b2ea8806ced 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -64,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); @@ -90,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/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 9bd057f8cd1..05db78ac40b 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -71,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); @@ -101,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 69cfa1a8073..93553b0bc29 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -75,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); @@ -104,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 fb821adaa03..85fd70a3520 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -185,11 +185,11 @@ void Reactor::updateConnected(bool updatePressure) { // Update the mass flow rate of connected flow devices double time = m_net->time(); - for (auto device : m_inlet) { - device->updateMassFlowRate(time); + for (size_t i = 0; i < m_outlet.size(); i++) { + m_outlet[i]->updateMassFlowRate(time); } - for (auto device : m_outlet) { - device->updateMassFlowRate(time); + for (size_t i = 0; i < m_inlet.size(); i++) { + m_inlet[i]->updateMassFlowRate(time); } } @@ -199,7 +199,6 @@ void Reactor::evalEqs(doublereal time, doublereal* y, double dmdt = 0.0; // dm/dt (gas phase) double* dYdt = ydot + 3; - evalFlowDevices(time); evalWalls(time); applySensitivity(params); m_thermo->restoreState(m_state); @@ -234,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(); } } @@ -269,18 +270,6 @@ void Reactor::evalWalls(double t) } } -void Reactor::evalFlowDevices(double t) -{ - for (size_t i = 0; i < m_outlet.size(); i++) { - m_outlet[i]->updateMassFlowRate(t); - m_mdot_out[i] = m_outlet[i]->massFlowRate(); - } - for (size_t i = 0; i < m_inlet.size(); i++) { - m_inlet[i]->updateMassFlowRate(t); - m_mdot_in[i] = m_inlet[i]->massFlowRate(); - } -} - 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);