diff --git a/include/cantera/clib/ctreactor.h b/include/cantera/clib/ctreactor.h index 8f550b8bdf..f2bcfaa2e8 100644 --- a/include/cantera/clib/ctreactor.h +++ b/include/cantera/clib/ctreactor.h @@ -51,7 +51,9 @@ extern "C" { CANTERA_CAPI int flowdev_new(const char* type); CANTERA_CAPI int flowdev_del(int i); CANTERA_CAPI int flowdev_install(int i, int n, int m); + //! @deprecated To be removed after Cantera 3.0; replaced by flowdev_setPrimary CANTERA_CAPI int flowdev_setMaster(int i, int n); + CANTERA_CAPI int flowdev_setPrimary(int i, int n); CANTERA_CAPI double flowdev_massFlowRate(int i); CANTERA_CAPI int flowdev_setMassFlowCoeff(int i, double v); CANTERA_CAPI int flowdev_setValveCoeff(int i, double v); @@ -62,8 +64,12 @@ extern "C" { CANTERA_CAPI int wall_new(const char* type); CANTERA_CAPI int wall_del(int i); CANTERA_CAPI int wall_install(int i, int n, int m); + //! @deprecated Only used by traditional MATLAB toolbox CANTERA_CAPI double wall_vdot(int i, double t); + CANTERA_CAPI double wall_expansionRate(int i); + //! @deprecated Only used by traditional MATLAB toolbox CANTERA_CAPI double wall_Q(int i, double t); + CANTERA_CAPI double wall_heatRate(int i); CANTERA_CAPI double wall_area(int i); CANTERA_CAPI int wall_setArea(int i, double v); CANTERA_CAPI int wall_setThermalResistance(int i, double rth); diff --git a/include/cantera/numerics/Func1.h b/include/cantera/numerics/Func1.h index 6866ef5de1..63fbbd3b2e 100644 --- a/include/cantera/numerics/Func1.h +++ b/include/cantera/numerics/Func1.h @@ -132,7 +132,7 @@ class Func1 return "functor"; } - //! Returns a string describing the type of the function + //! Returns a string with the class name of the functor //! @since New in Cantera 3.0. string typeName() const; diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 2a31362afe..bf1e0feb22 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -77,16 +77,43 @@ class FlowDevice return *m_out; } + //! Return current value of the pressure function. + /*! + * The mass flow rate [kg/s] is calculated given the pressure drop [Pa] and a + * coefficient set by a flow device specific function; unless a user-defined + * pressure function is set, this is the pressure difference across the device. + * The calculation of mass flow rate depends to the flow device. + * @since New in Cantera 3.0. + */ + double evalPressureFunction(); + //! Set a function of pressure that is used in determining the //! mass flow rate through the device. The evaluation of mass flow //! depends on the derived flow device class. virtual void setPressureFunction(Func1* f); + //! Return current value of the time function. + /*! + * The mass flow rate [kg/s] is calculated for a Flow device, and multiplied by a + * function of time, which returns 1.0 unless a user-defined function is provided. + * The calculation of mass flow rate depends on the flow device. + * @since New in Cantera 3.0. + */ + double evalTimeFunction(); + //! Set a function of time that is used in determining //! the mass flow rate through the device. The evaluation of mass flow //! depends on the derived flow device class. virtual void setTimeFunction(Func1* g); + //! Set current reactor network time + /*! + * @since New in Cantera 3.0. + */ + void setSimTime(double time) { + m_time = time; + } + protected: double m_mdot = Undef; @@ -99,6 +126,9 @@ class FlowDevice //! Coefficient set by derived classes; used by updateMassFlowRate double m_coeff = 1.0; + //! Current reactor network time + double m_time = 0.; + private: size_t m_nspin = 0; size_t m_nspout = 0; diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index d09472b852..326dbabe49 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -261,8 +261,8 @@ 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(). + //! Update the state information needed by connected reactors, flow devices, + //! and reactor walls. 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 diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h index 37e2ba44b3..4996551080 100644 --- a/include/cantera/zeroD/ReactorDelegator.h +++ b/include/cantera/zeroD/ReactorDelegator.h @@ -29,21 +29,21 @@ class ReactorAccessor virtual void setNEq(size_t n) = 0; //! Get the net rate of volume change (for example, from moving walls) [m^3/s] - virtual double vdot() const = 0; + virtual double expansionRate() const = 0; //! Set the net rate of volume change (for example, from moving walls) [m^3/s] - virtual void setVdot(double v) = 0; + virtual void setExpansionRate(double v) = 0; //! Get the net heat transfer rate (for example, through walls) into the //! reactor [W]. This value is initialized and calculated as part of //! Reactor::evalWalls(). - virtual double qdot() const = 0; + virtual double heatRate() const = 0; //! Set the net heat transfer rate (for example, through walls) into the //! reactor [W]. For a value set using this method to affect the calculations done //! by Reactor::eval, this method should be called in either a "replace" or "after" //! delegate for Reactor::evalWalls(). - virtual void setQdot(double q) = 0; + virtual void setHeatRate(double q) = 0; //! Set the state of the thermo object to correspond to the state of the reactor virtual void restoreThermoState() = 0; @@ -160,19 +160,19 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor R::m_nv = n; } - virtual double vdot() const override { + virtual double expansionRate() const override { return R::m_vdot; } - virtual void setVdot(double v) override { + virtual void setExpansionRate(double v) override { R::m_vdot = v; } - virtual double qdot() const override { + virtual double heatRate() const override { return R::m_Qdot; } - virtual void setQdot(double q) override { + virtual void setHeatRate(double q) override { R::m_Qdot = q; } diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 3a819df3ee..06d9399b94 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -51,6 +51,14 @@ class ReactorNet : public FuncEval //! state as the initial condition. void setInitialTime(double time); + //! Get the initial value of the independent variable (typically time). + /*! + * @since New in Cantera 3.0. + */ + double getInitialTime() const { + return m_initial_time; + } + //! Get the maximum integrator step. double maxTimeStep() { return m_maxstep; @@ -316,6 +324,9 @@ class ReactorNet : public FuncEval //! on the type of reactors in the network. double m_time = 0.0; + //! The initial value of the independent variable in the system. + double m_initial_time = 0.0; + bool m_init = false; bool m_integrator_init = false; //!< True if integrator initialization is current size_t m_nv = 0; diff --git a/include/cantera/zeroD/Wall.h b/include/cantera/zeroD/Wall.h index 17fca0d833..50ee472a1b 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -40,8 +40,22 @@ class WallBase /*! * This method is called by Reactor::evalWalls(). Base class method * does nothing (that is, constant volume), but may be overloaded. + * @deprecated To be removed after Cantera 3.0; replaceable by expansionRate. */ virtual double vdot(double t) { + warn_deprecated("WallBase::vdot", + "To be removed after Cantera 3.0; replaceable by 'expansionRate'."); + return 0.0; + } + + //! Rate of volume change (m^3/s) for the adjacent reactors at current reactor + //! network time. + /*! + * This method is called by Reactor::evalWalls(). Base class method + * does nothing (that is, constant volume), but may be overloaded. + * @since New in Cantera 3.0. + */ + virtual double expansionRate() { return 0.0; } @@ -49,8 +63,21 @@ class WallBase /*! * This method is called by Reactor::evalWalls(). Base class method * does nothing (that is, an adiabatic wall), but may be overloaded. + * @deprecated To be removed after Cantera 3.0; replaceable by heatRate. */ virtual double Q(double t) { + warn_deprecated("WallBase::Q", + "To be removed after Cantera 3.0; replaceable by 'heatRate'."); + return 0.0; + } + + //! Heat flow rate through the wall (W) at current reactor network time. + /*! + * This method is called by Reactor::evalWalls(). Base class method + * does nothing (that is, an adiabatic wall), but may be overloaded. + * @since New in Cantera 3.0. + */ + virtual double heatRate() { return 0.0; } @@ -83,12 +110,23 @@ class WallBase return *m_right; } + //! Set current reactor network time + /*! + * @since New in Cantera 3.0. + */ + void setSimTime(double time) { + m_time = time; + } + protected: ReactorBase* m_left = nullptr; ReactorBase* m_right = nullptr; std::vector m_surf; + //! current reactor network time + double m_time = 0.0; + double m_area = 1.0; }; @@ -108,6 +146,10 @@ class Wall : public WallBase return "Wall"; } + //! Wall velocity \f$ v(t) \f$ at current reactor network time. + //! @since New in Cantera 3.0. + double velocity() const; + //! Set the wall velocity to a specified function of time, \f$ v(t) \f$. void setVelocity(Func1* f=0) { if (f) { @@ -125,9 +167,29 @@ class Wall : public WallBase * area, and *F(t)* is a specified function of time. Positive values for * `vdot` correspond to increases in the volume of reactor on left, and * decreases in the volume of the reactor on the right. + * @deprecated Still used by traditional MATLAB toolbox; replaceable by + * expansionRate. */ virtual double vdot(double t); + //! Rate of volume change (m^3/s) for the adjacent reactors. + /*! + * The volume rate of change is given by + * \f[ + * \dot V = K A (P_{left} - P_{right}) + F(t) + * \f] + * where *K* is the specified expansion rate coefficient, *A* is the wall area, + * and and *F(t)* is a specified function evaluated at the current network time. + * Positive values for `expansionRate` correspond to increases in the volume of + * reactor on left, and decreases in the volume of the reactor on the right. + * @since New in Cantera 3.0. + */ + virtual double expansionRate(); + + //! Heat flux function \f$ q_0(t) \f$ evaluated at current reactor network time. + //! @since New in Cantera 3.0. + double heatFlux() const; + //! Specify the heat flux function \f$ q_0(t) \f$. void setHeatFlux(Func1* q) { m_qf = q; @@ -142,9 +204,23 @@ class Wall : public WallBase * where *h* is the heat transfer coefficient, *A* is the wall area, and * *G(t)* is a specified function of time. Positive values denote a flux * from left to right. + * @deprecated Still used by traditional MATLAB toolbox; replaceable by heatRate. */ virtual double Q(double t); + //! Heat flow rate through the wall (W). + /*! + * The heat flux is given by + * \f[ + * Q = h A (T_{left} - T_{right}) + A G(t) + * \f] + * where *h* is the heat transfer coefficient, *A* is the wall area, and + * *G(t)* is a specified function of time evaluated at the current network + * time. Positive values denote a flux from left to right. + * @since New in Cantera 3.0. + */ + virtual double heatRate(); + void setThermalResistance(double Rth) { m_rrth = 1.0/Rth; } diff --git a/include/cantera/zeroD/flowControllers.h b/include/cantera/zeroD/flowControllers.h index a0763a144c..840489a3a6 100644 --- a/include/cantera/zeroD/flowControllers.h +++ b/include/cantera/zeroD/flowControllers.h @@ -57,7 +57,7 @@ class MassFlowController : public FlowDevice /** * A class for flow controllers where the flow rate is equal to the flow rate - * of a "master" mass flow controller plus a correction proportional to the + * of a primary mass flow controller plus a correction proportional to the * pressure difference between the inlet and outlet. */ class PressureController : public FlowDevice @@ -69,15 +69,21 @@ class PressureController : public FlowDevice return "PressureController"; } - virtual bool ready() { - return FlowDevice::ready() && m_master != 0; + return FlowDevice::ready() && m_primary != 0; } - void setMaster(FlowDevice* master) { - m_master = master; + //! Set the primary mass flow controller. + /*! + * @since New in Cantera 3.0. + */ + void setPrimary(FlowDevice* primary) { + m_primary = primary; } + //! @deprecated To be removed after Cantera 3.0; replaced by setPrimary. + void setMaster(FlowDevice* master); + virtual void setTimeFunction(Func1* g) { throw NotImplementedError("PressureController::setTimeFunction"); } @@ -86,11 +92,11 @@ class PressureController : public FlowDevice //! rate /*! * *c* has units of kg/s/Pa. The mass flow rate is computed as: - * \f[\dot{m} = \dot{m}_{master} + c f(\Delta P) \f] + * \f[\dot{m} = \dot{m}_{primary} + c f(\Delta P) \f] * where *f* is a functions of pressure drop that is set by * `setPressureFunction`. If no functions is specified, the mass flow * rate defaults to: - * \f[\dot{m} = \dot{m}_{master} + c \Delta P \f] + * \f[\dot{m} = \dot{m}_{primary} + c \Delta P \f] */ void setPressureCoeff(double c) { m_coeff = c; @@ -104,7 +110,7 @@ class PressureController : public FlowDevice virtual void updateMassFlowRate(double time); protected: - FlowDevice* m_master = nullptr; + FlowDevice* m_primary = nullptr; }; //! Supply a mass flow rate that is a function of the pressure drop across the diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index bfc0705163..6baac1f1a4 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -94,8 +94,10 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": void setCoverages(int, double*) void setCoverages(int, Composition&) except +translate_exception void syncCoverages(int) - double vdot(double) - double Q(double) + double expansionRate() except +translate_exception + double vdot(double) except +translate_exception + double heatRate() except +translate_exception + double Q(double) except +translate_exception void addSensitivityReaction(int, size_t) except +translate_exception size_t nSensParams(int) @@ -108,7 +110,9 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": double getHeatTransferCoeff() void setEmissivity(double) except +translate_exception double getEmissivity() + double velocity() void setVelocity(CxxFunc1*) + double heatFlux() void setHeatFlux(CxxFunc1*) # reactor surface @@ -132,7 +136,9 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": double massFlowRate() except +translate_exception double massFlowRate(double) except +translate_exception cbool install(CxxReactorBase&, CxxReactorBase&) except +translate_exception + double evalPressureFunction() except +translate_exception void setPressureFunction(CxxFunc1*) except +translate_exception + double evalTimeFunction() except +translate_exception void setTimeFunction(CxxFunc1*) except +translate_exception cdef cppclass CxxMassFlowController "Cantera::MassFlowController" (CxxFlowDevice): @@ -150,7 +156,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": CxxPressureController() double getPressureCoeff() void setPressureCoeff(double) - void setMaster(CxxFlowDevice*) + void setPrimary(CxxFlowDevice*) # reactor net @@ -164,6 +170,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": double time() except +translate_exception double distance() except +translate_exception void setInitialTime(double) + double getInitialTime() void setTolerances(double, double) double rtol() double atol() @@ -199,10 +206,10 @@ cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera": cdef cppclass CxxReactorAccessor "Cantera::ReactorAccessor": CxxReactorAccessor() void setNEq(size_t) - double vdot() - void setVdot(double) - double qdot() - void setQdot(double) + double expansionRate() + void setExpansionRate(double) + double heatRate() + void setHeatRate(double) void restoreThermoState() except +translate_exception void restoreSurfaceState(size_t) except +translate_exception diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 785f7abd62..e6bcb35c35 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -595,8 +595,8 @@ cdef class ExtensibleReactor(Reactor): `n_vars`. ``eval_walls(self, t : double) -> None`` - Responsible for calculating the net rate of volume change `vdot` - and the net rate of heat transfer `qdot` caused by walls connected + Responsible for calculating the net rate of volume change `expansion_rate` + and the net rate of heat transfer `heat_rate` caused by walls connected to this reactor. ``eval_surfaces(LHS : double[:], RHS : double[:], sdot : double[:]) -> None`` @@ -654,20 +654,66 @@ cdef class ExtensibleReactor(Reactor): property vdot: """ Get/Set the net rate of volume change (for example, from moving walls) [m^3/s] + + .. deprecated:: 3.0 + + To be removed in Cantera 3.0; renamed to `expansion_rate`. """ def __get__(self): - return self.accessor.vdot() + warnings.warn( + "ExtensibleReactor.vdot: To be removed in Cantera 3.0; " + "renamed to 'expansion_rate'.", DeprecationWarning) + return self.accessor.expansionRate() def __set__(self, vdot): - self.accessor.setVdot(vdot) + warnings.warn( + "ExtensibleReactor.vdot: To be removed in Cantera 3.0; " + "renamed to 'expansion_rate'.", DeprecationWarning) + self.accessor.setExpansionRate(vdot) + + @property + def expansion_rate(self): + """ + Get/Set the net rate of volume change (for example, from moving walls) [m^3/s] + + .. versionadded:: 3.0 + """ + return self.accessor.expansionRate() + + @expansion_rate.setter + def expansion_rate(self, vdot): + self.accessor.setExpansionRate(vdot) property qdot: """ Get/Set the net heat transfer rate (for example, through walls) [W] + + .. deprecated:: 3.0 + + To be removed in Cantera 3.0; renamed to `heat_rate`. """ def __get__(self): - return self.accessor.qdot() + warnings.warn( + "ExtensibleReactor.qdot: To be removed in Cantera 3.0; " + "renamed to 'heat_rate'.", DeprecationWarning) + return self.accessor.heatRate() def __set__(self, qdot): - self.accessor.setQdot(qdot) + warnings.warn( + "ExtensibleReactor.qdot: To be removed in Cantera 3.0; " + "renamed to 'heat_rate'.", DeprecationWarning) + self.accessor.setHeatRate(qdot) + + @property + def heat_rate(self): + """ + Get/Set the net heat transfer rate (for example, through walls) [W] + + .. versionadded:: 3.0 + """ + return self.accessor.heatRate() + + @heat_rate.setter + def heat_rate(self, qdot): + self.accessor.setHeatRate(qdot) def restore_thermo_state(self): """ @@ -883,9 +929,9 @@ cdef class WallBase: if U is not None: self.heat_transfer_coeff = U if Q is not None: - self.set_heat_flux(Q) + self.heat_flux = Q if velocity is not None: - self.set_velocity(velocity) + self.velocity = velocity def _install(self, ReactorBase left, ReactorBase right): """ @@ -911,19 +957,49 @@ cdef class WallBase: def __set__(self, double value): self.wall.setArea(value) + @property + def expansion_rate(self): + """ + Get the rate of volumetric change [m^3/s] associated with the wall at the + current reactor network time. A positive value corresponds to the left-hand + reactor volume increasing, and the right-hand reactor volume decreasing. + + .. versionadded:: 3.0 + """ + return self.wall.expansionRate() + def vdot(self, double t): """ The rate of volumetric change [m^3/s] associated with the wall at time ``t``. A positive value corresponds to the left-hand reactor volume increasing, and the right-hand reactor volume decreasing. + + .. deprecated:: 3.0 + + To be removed after Cantera 3.0; replaceable by ``expansion_rate``. """ return self.wall.vdot(t) + @property + def heat_rate(self): + """ + Get the total heat flux [W] through the wall at the current reactor network + time. A positive value corresponds to heat flowing from the left-hand reactor + to the right-hand one. + + .. versionadded:: 3.0 + """ + return self.wall.heatRate() + def qdot(self, double t): """ Total heat flux [W] through the wall at time ``t``. A positive value corresponds to heat flowing from the left-hand reactor to the right-hand one. + + .. deprecated:: 3.0 + + To be removed after Cantera 3.0; replaceable by ``heat_rate``. """ return self.wall.Q(t) @@ -981,11 +1057,18 @@ cdef class Wall(WallBase): def __set__(self, double value): ((self.wall)).setEmissivity(value) - def set_velocity(self, v): + @property + def velocity(self): """ - The wall velocity [m/s]. May be either a constant or an arbitrary + The wall velocity [m/s]. May be either set to a constant or an arbitrary function of time. See `Func1`. + + .. versionadded:: 3.0 """ + return ((self.wall)).velocity() + + @velocity.setter + def velocity(self, v): cdef Func1 f if isinstance(v, Func1): f = v @@ -995,11 +1078,32 @@ cdef class Wall(WallBase): self._velocity_func = f ((self.wall)).setVelocity(f.func) - def set_heat_flux(self, q): + def set_velocity(self, v): """ - Heat flux [W/m^2] across the wall. May be either a constant or + The wall velocity [m/s]. May be either a constant or an arbitrary + function of time. See `Func1`. + + .. deprecated:: 3.0 + + Replaced by the ``velocity`` property. + """ + warnings.warn( + "Wall.set_velocity: To be removed after Cantera 3.0; replaced by property " + "'velocity'.", DeprecationWarning) + self.velocity = v + + @property + def heat_flux(self): + """ + Heat flux [W/m^2] across the wall. May be either set to a constant or an arbitrary function of time. See `Func1`. + + .. versionadded:: 3.0 """ + return ((self.wall)).heatFlux() + + @heat_flux.setter + def heat_flux(self, q): cdef Func1 f if isinstance(q, Func1): f = q @@ -1009,6 +1113,20 @@ cdef class Wall(WallBase): self._heat_flux_func = f (self.wall).setHeatFlux(f.func) + def set_heat_flux(self, q): + """ + Heat flux [W/m^2] across the wall. May be either a constant or + an arbitrary function of time. See `Func1`. + + .. deprecated:: 3.0 + + Replaced by the ``heat_flux`` property. + """ + warnings.warn( + "Wall.set_heat_flux: To be removed after Cantera 3.0; replaced by property " + "'heat_flux'.", DeprecationWarning) + self.heat_flux = q + cdef class FlowDevice: """ @@ -1064,6 +1182,34 @@ cdef class FlowDevice: def __get__(self): return self.dev.massFlowRate() + @property + def pressure_function(self): + r""" + The relationship between mass flow rate and the pressure drop across a flow + device. The mass flow rate [kg/s] is calculated given the pressure drop [Pa] and + a coefficient set by a flow device specific function. Unless a user-defined + pressure function is provided, the function returns the pressure difference + across the device. The calculation of mass flow rate depends on the flow device. + + >>> f = FlowDevice(res1, reactor1) + >>> f.pressure_function = lambda dP: dP**2 + + where `FlowDevice` is either a `Valve` or `PressureController` object. + + .. versionadded:: 3.0 + """ + return self.dev.evalPressureFunction() + + @pressure_function.setter + def pressure_function(self, k): + cdef Func1 f + if isinstance(k, Func1): + f = k + else: + f = Func1(k) + self._rate_func = f + self.dev.setPressureFunction(f.func) + def set_pressure_function(self, k): r""" Set the relationship between mass flow rate and the pressure drop across a @@ -1075,14 +1221,41 @@ cdef class FlowDevice: >>> F.set_pressure_function(lambda dP: dP**2) where FlowDevice is either a Valve or PressureController object. + + .. deprecated:: 3.0 + To be removed after Cantera 3.0. Use property ``pressure_function`` instead. """ - cdef Func1 f + warnings.warn( + "FlowDevice.set_pressure_function: To be removed after Cantera 3.0; " + "replaced by 'pressure_function'.", DeprecationWarning) + self.pressure_function = k + + @property + def time_function(self): + r""" + The time dependence of a flow device. The mass flow rate [kg/s] is calculated + for a Flow device, and multiplied by a function of time, which returns 1.0 + unless a user-defined function is provided. The calculation of mass flow rate + depends on the flow device. + + >>> f = FlowDevice(res1, reactor1) + >>> f.time_function = lambda t: exp(-10 * (t - 0.5)**2) + + where `FlowDevice` is either a `Valve` or `MassFlowController` object. + + .. versionadded:: 3.0 + """ + return self.dev.evalTimeFunction() + + @time_function.setter + def time_function(self, k): + cdef Func1 g if isinstance(k, Func1): - f = k + g = k else: - f = Func1(k) - self._rate_func = f - self.dev.setPressureFunction(f.func) + g = Func1(k) + self._time_func = g + self.dev.setTimeFunction(g.func) def set_time_function(self, k): r""" @@ -1094,14 +1267,14 @@ cdef class FlowDevice: >>> F.set_time_function(lambda t: exp(-10 * (t - 0.5)**2)) where FlowDevice is either a Valve or MassFlowController object. + + .. deprecated:: 3.0 + To be removed after Cantera 3.0. Use property ``time_function`` instead. """ - cdef Func1 g - if isinstance(k, Func1): - g = k - else: - g = Func1(k) - self._time_func = g - self.dev.setTimeFunction(g.func) + warnings.warn( + "FlowDevice.set_time_function: To be removed after Cantera 3.0; " + "replaced by 'time_function'.", DeprecationWarning) + self.time_function = k cdef class MassFlowController(FlowDevice): @@ -1114,10 +1287,10 @@ cdef class MassFlowController(FlowDevice): where :math:`\dot m_0` is a constant value and :math:`g(t)` is a function of time. Both :math:`\dot m_0` and :math:`g(t)` can be set individually by - the property `mass_flow_coeff` and the method `set_time_function`, - respectively. The property `mass_flow_rate` combines the former - into a single interface. Note that if :math:`\dot m_0*g(t) < 0`, the mass flow - rate will be set to zero, since reversal of the flow direction is not allowed. + properties `mass_flow_coeff` and `time_function`, respectively. The property + `mass_flow_rate` combines the former into a single interface. Note that if + :math:`\dot m_0*g(t) < 0`, the mass flow rate will be set to zero, since + reversal of the flow direction is not allowed. Unlike a real mass flow controller, a MassFlowController object will maintain the flow even if the downstream pressure is greater than the @@ -1135,7 +1308,7 @@ cdef class MassFlowController(FlowDevice): property mass_flow_coeff: r"""Set the mass flow rate [kg/s] through the mass flow controller as a constant, which may be modified by a function of time, see - `set_time_function`. + `time_function`. >>> mfc = MassFlowController(res1, reactor1) >>> mfc.mass_flow_coeff = 1e-4 # Set the flow rate to a constant @@ -1153,7 +1326,7 @@ cdef class MassFlowController(FlowDevice): 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. + the property `mass_flow_coeff` or updates the `time_function` property. >>> mfc.mass_flow_rate = 0.3 >>> mfc.mass_flow_rate = lambda t: 2.5 * exp(-10 * (t - 0.5)**2) @@ -1166,7 +1339,7 @@ cdef class MassFlowController(FlowDevice): (self.dev).setMassFlowRate(m) else: self.mass_flow_coeff = 1. - self.set_time_function(m) + self.time_function = m cdef class Valve(FlowDevice): @@ -1185,12 +1358,12 @@ cdef class Valve(FlowDevice): where :math:`f` is the arbitrary function that multiplies :math:`K_v` given a single argument, the pressure differential. Further, a valve opening function - :math:`g` may be specified using the method `set_time_function`, such that + :math:`g` may be specified using the `time_function` property, such that .. math:: \dot m = K_v*g(t)*f(P_1 - P_2) See the documentation for the `valve_coeff` property as well as the - `set_pressure_function` and `set_time_function` methods for examples. Note that + `pressure_function` and `time_function` properties for examples. Note that it is never possible for the flow to reverse and go from the downstream to the upstream reactor/reservoir through a line containing a `Valve` object. @@ -1208,15 +1381,15 @@ cdef class Valve(FlowDevice): self.valve_coeff = K else: self.valve_coeff = 1. - self.set_pressure_function(K) + self.pressure_function = K property valve_coeff: r"""Set valve coefficient, that is, the proportionality constant between mass flow rate and pressure drop [kg/s/Pa]. - >>> V = Valve(res1, reactor1) - >>> V.valve_coeff = 1e-4 # Set the value of K to a constant - >>> V.valve_coeff # Get the value of K + >>> v = Valve(res1, reactor1) + >>> v.valve_coeff = 1e-4 # Set the value of K to a constant + >>> v.valve_coeff # Get the value of K """ def __get__(self): return (self.dev).getValveCoeff() @@ -1227,33 +1400,39 @@ cdef class Valve(FlowDevice): cdef class PressureController(FlowDevice): r""" A PressureController is designed to be used in conjunction with another - 'master' flow controller, typically a `MassFlowController`. The master + primary flow controller, typically a `MassFlowController`. The primary flow controller is installed on the inlet of the reactor, and the corresponding `PressureController` is installed on the outlet of the - reactor. The `PressureController` mass flow rate is equal to the master + reactor. The `PressureController` mass flow rate is equal to the primary mass flow rate, plus a small correction dependent on the pressure difference: - .. math:: \dot m = \dot m_{\rm master} + K_v(P_1 - P_2). + .. math:: \dot m = \dot m_{\rm primary} + K_v(P_1 - P_2). As an alternative, an arbitrary function of pressure differential can be - specified using the method `set_pressure_function`, such that + specified using the `pressure_function` property, such that - .. math:: \dot m = \dot m_{\rm master} + K_v*f(P_1 - P_2) + .. math:: \dot m = \dot m_{\rm primary} + K_v*f(P_1 - P_2) where :math:`f` is the arbitrary function of a single argument. """ flowdevice_type = "PressureController" - def __init__(self, upstream, downstream, *, name=None, master=None, K=1.): + def __init__(self, upstream, downstream, *, + name=None, primary=None, K=1., **kwargs): + if "master" in kwargs: + warnings.warn( + "PressureController: The 'master' keyword argument is deprecated; " + "use 'primary' instead.", DeprecationWarning) + primary = kwargs["master"] super().__init__(upstream, downstream, name=name) - if master is not None: - self.set_master(master) + if primary is not None: + self.primary = primary if isinstance(K, _numbers.Real): self.pressure_coeff = K else: self.pressure_coeff = 1. - self.set_pressure_function(K) + self.pressure_function = K property pressure_coeff: """ @@ -1265,12 +1444,32 @@ cdef class PressureController(FlowDevice): def __set__(self, double value): (self.dev).setPressureCoeff(value) + @property + def primary(self): + """ + Primary `FlowDevice` used to compute this device's mass flow rate. + + .. versionadded:: 3.0 + """ + raise NotImplementedError("PressureController.primary") + + @primary.setter + def primary(self, FlowDevice d): + (self.dev).setPrimary(d.dev) + def set_master(self, FlowDevice d): """ Set the "master" `FlowDevice` used to compute this device's mass flow rate. + + .. deprecated:: 3.0 + + To be removed after Cantera 3.0; replaced by property ``primary``. """ - (self.dev).setMaster(d.dev) + warnings.warn( + "PressureController.set_master: To be removed after Cantera 3.0; " + "replaced by 'primary'.", DeprecationWarning) + self.primary = d cdef class ReactorNet: @@ -1348,12 +1547,32 @@ cdef class ReactorNet: """ return self.net.distance() + @property + def initial_time(self): + """ + The initial time of the integrator. When set, integration is restarted from this + time using the current state as the initial condition. Default: 0.0 s. + + .. versionadded:: 3.0 + """ + return self.net.getInitialTime() + + @initial_time.setter + def initial_time(self, double t): + self.net.setInitialTime(t) + def set_initial_time(self, double t): """ Set the initial time. Restarts integration from this time using the current state as the initial condition. Default: 0.0 s. + + .. deprecated:: 3.0 + To be removed after Cantera 3.0. Use property ``initial_time`` instead. """ - self.net.setInitialTime(t) + warnings.warn( + "ReactorNet.set_initial_time: To be removed after Cantera 3.0. " + "Use property 'initial_time' instead.", DeprecationWarning) + self.initial_time = t property max_time_step: """ diff --git a/interfaces/matlab_experimental/Reactor/FlowDevice.m b/interfaces/matlab_experimental/Reactor/FlowDevice.m index 62d47b65bd..ca89e1998d 100644 --- a/interfaces/matlab_experimental/Reactor/FlowDevice.m +++ b/interfaces/matlab_experimental/Reactor/FlowDevice.m @@ -134,11 +134,11 @@ function install(f, upstream, downstream) end - function setMaster(f, d) - % Set the Master flow device used to compute this device's + function setPrimary(f, d) + % Set the primary flow device used to compute this device's % mass flow rate. :: % - % >> f.setMaster(d) + % >> f.setPrimary(d) % % :param f: % Instance of class :mat:class:`MassFlowController`. @@ -146,9 +146,9 @@ function setMaster(f, d) % Instance of class :mat:class:`Func`. if strcmp(f.type, 'PressureController') - k = ctFunc('flowdev_setMaster', f.id, d); + k = ctFunc('flowdev_setPrimary', f.id, d); else - error('Master flow device can only be set for pressure controllers.'); + error('Primary flow device can only be set for pressure controllers.'); end end diff --git a/interfaces/matlab_experimental/Reactor/Wall.m b/interfaces/matlab_experimental/Reactor/Wall.m index ce78f11565..ecaa3e47db 100644 --- a/interfaces/matlab_experimental/Reactor/Wall.m +++ b/interfaces/matlab_experimental/Reactor/Wall.m @@ -60,6 +60,8 @@ properties (SetAccess = public) area % Area of the wall in m^2. + heatRate % Total heat transfer rate through the wall at current time step in W. + expansionRate % Rate of volumetric change at current time step in m^3/s. thermalResistance % Thermal resistance in K*m^2/W. heatTransferCoeff % Heat transfer coefficient in W/(m^2-K). emissivity % Non-dimensional emissivity. @@ -129,15 +131,12 @@ function delete(w) a = ctFunc('wall_area', w.id); end - function q = qdot(w, t) - % Total heat transfer rate through a wall at a given time t. - - q = ctFunc('wall_Q', w.id, t); + function q = get.heatRate(w) + q = ctFunc('wall_heatRate', w.id); end - function v = vdot(w, t) - % Rate of volumetric change at a given time t. - v = ctFunc('wall_vdot', w.id, t); + function v = get.expansionRate(w) + v = ctFunc('wall_expansionRate', w.id); end %% ReactorNet set methods diff --git a/samples/python/reactors/combustor.py b/samples/python/reactors/combustor.py index 08e12d5595..89df4c5bd4 100644 --- a/samples/python/reactors/combustor.py +++ b/samples/python/reactors/combustor.py @@ -11,7 +11,7 @@ enclosing scope. Also shows the use of a PressureController to create a constant pressure reactor with a fixed volume. -Requires: cantera >= 2.5.0, matplotlib >= 2.0 +Requires: cantera >= 3.0, matplotlib >= 2.0 Keywords: combustion, reactor network, well-stirred reactor, plotting """ @@ -54,11 +54,11 @@ def mdot(t): inlet_mfc = ct.MassFlowController(inlet, combustor, mdot=mdot) -# A PressureController has a baseline mass flow rate matching the 'master' +# A PressureController has a baseline mass flow rate matching the 'primary' # MassFlowController, with an additional pressure-dependent term. By explicitly # including the upstream mass flow rate, the pressure is kept constant without # needing to use a large value for 'K', which can introduce undesired stiffness. -outlet_mfc = ct.PressureController(combustor, exhaust, master=inlet_mfc, K=0.01) +outlet_mfc = ct.PressureController(combustor, exhaust, primary=inlet_mfc, K=0.01) # the simulation only contains one reactor sim = ct.ReactorNet([combustor]) @@ -69,7 +69,7 @@ def mdot(t): residence_time = 0.1 # starting residence time while combustor.T > 500: - sim.set_initial_time(0.0) # reset the integrator + sim.initial_time = 0.0 # reset the integrator sim.advance_to_steady_state() print('tres = {:.2e}; T = {:.1f}'.format(residence_time, combustor.T)) states.append(combustor.thermo.state, tres=residence_time) diff --git a/samples/python/reactors/custom2.py b/samples/python/reactors/custom2.py index c418df53bb..9c519f3274 100644 --- a/samples/python/reactors/custom2.py +++ b/samples/python/reactors/custom2.py @@ -14,7 +14,7 @@ determined by integrating the equation of motion. This requires adding a new variable to the reactor's state vector which represents the wall velocity. -Requires: cantera >= 2.6.0, matplotlib >= 2.0 +Requires: cantera >= 3.0, matplotlib >= 2.0 Keywords: combustion, reactor network, user-defined model, plotting """ @@ -47,7 +47,7 @@ def after_update_state(self, y): # This method is used to set the state of the Reactor and Wall objects # based on the new values for the state vector provided by the ODE solver self.v_wall = y[self.i_wall] - self.walls[0].set_velocity(self.v_wall) + self.walls[0].velocity = self.v_wall def after_eval(self, t, LHS, RHS): # Calculate the time derivative for the additional equation diff --git a/samples/python/reactors/ic_engine.py b/samples/python/reactors/ic_engine.py index caf5c93976..e7978b6409 100644 --- a/samples/python/reactors/ic_engine.py +++ b/samples/python/reactors/ic_engine.py @@ -1,4 +1,3 @@ -# -*- coding: utf-8 -*- """ Simulation of a (gaseous) Diesel-type internal combustion engine. @@ -6,7 +5,7 @@ center. Note that this example uses numerous simplifying assumptions and thus serves for illustration purposes only. -Requires: cantera >= 2.5.0, scipy >= 0.19, matplotlib >= 2.0 +Requires: cantera >= 3.0, scipy >= 0.19, matplotlib >= 2.0 Keywords: combustion, thermodynamics, internal combustion engine, thermodynamic cycle, reactor network, plotting, pollutant formation """ @@ -104,13 +103,15 @@ def piston_speed(t): # define inlet state gas.TPX = T_inlet, p_inlet, comp_inlet +# Note: The previous line is technically not needed as the state of the gas object is +# already set correctly; change if inlet state is different from the reactor state. inlet = ct.Reservoir(gas) # inlet valve inlet_valve = ct.Valve(inlet, cyl) inlet_delta = np.mod(inlet_close - inlet_open, 4 * np.pi) inlet_valve.valve_coeff = inlet_valve_coeff -inlet_valve.set_time_function( +inlet_valve.time_function = ( lambda t: np.mod(crank_angle(t) - inlet_open, 4 * np.pi) < inlet_delta) # define injector state (gaseous!) @@ -122,7 +123,7 @@ def piston_speed(t): injector_delta = np.mod(injector_close - injector_open, 4 * np.pi) injector_t_open = (injector_close - injector_open) / 2. / np.pi / f injector_mfc.mass_flow_coeff = injector_mass / injector_t_open -injector_mfc.set_time_function( +injector_mfc.time_function = ( lambda t: np.mod(crank_angle(t) - injector_open, 4 * np.pi) < injector_delta) # define outlet pressure (temperature and composition don't matter) @@ -133,7 +134,7 @@ def piston_speed(t): outlet_valve = ct.Valve(cyl, outlet) outlet_delta = np.mod(outlet_close - outlet_open, 4 * np.pi) outlet_valve.valve_coeff = outlet_valve_coeff -outlet_valve.set_time_function( +outlet_valve.time_function = ( lambda t: np.mod(crank_angle(t) - outlet_open, 4 * np.pi) < outlet_delta) # define ambient pressure (temperature and composition don't matter) @@ -143,7 +144,7 @@ def piston_speed(t): # piston is modeled as a moving wall piston = ct.Wall(ambient_air, cyl) piston.area = A_piston -piston.set_velocity(piston_speed) +piston.velocity = piston_speed # create a reactor network containing the cylinder and limit advance step sim = ct.ReactorNet([cyl]) diff --git a/samples/python/reactors/pfr.py b/samples/python/reactors/pfr.py index abc8f3fb20..84ea319a28 100644 --- a/samples/python/reactors/pfr.py +++ b/samples/python/reactors/pfr.py @@ -1,10 +1,9 @@ -# -*- coding: utf-8 -*- """ This example solves a plug-flow reactor problem of hydrogen-oxygen combustion. The PFR is computed by two approaches: The simulation of a Lagrangian fluid particle, and the simulation of a chain of reactors. -Requires: cantera >= 2.5.0, matplotlib >= 2.0 +Requires: cantera >= 3.0, matplotlib >= 2.0 Keywords: combustion, reactor network, plug flow reactor """ @@ -109,7 +108,7 @@ # We need an outlet to the downstream reservoir. This will determine the # pressure in the reactor. The value of K will only affect the transient # pressure difference. -v = ct.PressureController(r2, downstream, master=m, K=1e-5) +v = ct.PressureController(r2, downstream, primary=m, K=1e-5) sim2 = ct.ReactorNet([r2]) diff --git a/samples/python/reactors/surf_pfr_chain.py b/samples/python/reactors/surf_pfr_chain.py index 3c8f803967..85e40e766f 100644 --- a/samples/python/reactors/surf_pfr_chain.py +++ b/samples/python/reactors/surf_pfr_chain.py @@ -5,7 +5,7 @@ DAE system, the PFR is approximated as a chain of successive WSRs. See surf_pfr.py for a more advanced implementation that solves the DAE system directly. -Requires: cantera >= 2.5.0 +Requires: cantera >= 3.0 Keywords: catalysis, reactor network, surface chemistry, plug flow reactor, packed bed reactor """ @@ -95,7 +95,7 @@ # We need an outlet to the downstream reservoir. This will determine the # pressure in the reactor. The value of K will only affect the transient # pressure difference. -v = ct.PressureController(r, downstream, master=m, K=1e-5) +v = ct.PressureController(r, downstream, primary=m, K=1e-5) sim = ct.ReactorNet([r]) diff --git a/src/clib/ctreactor.cpp b/src/clib/ctreactor.cpp index f9ccc56afa..1ef59ef9a7 100644 --- a/src/clib/ctreactor.cpp +++ b/src/clib/ctreactor.cpp @@ -394,6 +394,17 @@ extern "C" { } } + int flowdev_setPrimary(int i, int n) + { + try { + FlowDeviceCabinet::get(i).setPrimary( + &FlowDeviceCabinet::item(n)); + return 0; + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + double flowdev_massFlowRate(int i) { try { @@ -494,6 +505,15 @@ extern "C" { } } + double wall_expansionRate(int i) + { + try { + return WallCabinet::item(i).expansionRate(); + } catch (...) { + return handleAllExceptions(DERR, DERR); + } + } + double wall_Q(int i, double t) { try { @@ -503,6 +523,15 @@ extern "C" { } } + double wall_heatRate(int i) + { + try { + return WallCabinet::item(i).heatRate(); + } catch (...) { + return handleAllExceptions(DERR, DERR); + } + } + double wall_area(int i) { try { diff --git a/src/zeroD/FlowDevice.cpp b/src/zeroD/FlowDevice.cpp index 2e32921120..2384e860bd 100644 --- a/src/zeroD/FlowDevice.cpp +++ b/src/zeroD/FlowDevice.cpp @@ -47,11 +47,28 @@ void FlowDevice::setPressureFunction(Func1* f) m_pfunc = f; } +double FlowDevice::evalPressureFunction() +{ + double delta_P = in().pressure() - out().pressure(); + if (m_pfunc) { + return m_pfunc->eval(delta_P); + } + return delta_P; +} + void FlowDevice::setTimeFunction(Func1* g) { m_tfunc = g; } +double FlowDevice::evalTimeFunction() +{ + if (m_tfunc) { + return m_tfunc->eval(m_time); + } + return 1.; +} + double FlowDevice::outletSpeciesMassFlowRate(size_t k) { if (k >= m_nspout) { diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 762b08ffc6..014d9627d2 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -196,11 +196,16 @@ void Reactor::updateConnected(bool updatePressure) { time = (timeIsIndependent()) ? m_net->time() : m_net->distance(); } for (size_t i = 0; i < m_outlet.size(); i++) { + m_outlet[i]->setSimTime(time); m_outlet[i]->updateMassFlowRate(time); } for (size_t i = 0; i < m_inlet.size(); i++) { + m_inlet[i]->setSimTime(time); m_inlet[i]->updateMassFlowRate(time); } + for (size_t i = 0; i < m_wall.size(); i++) { + m_wall[i]->setSimTime(time); + } } void Reactor::eval(double time, double* LHS, double* RHS) @@ -269,12 +274,13 @@ void Reactor::eval(double time, double* LHS, double* RHS) void Reactor::evalWalls(double t) { + // time is currently unused m_vdot = 0.0; m_Qdot = 0.0; for (size_t i = 0; i < m_wall.size(); i++) { int f = 2 * m_lr[i] - 1; - m_vdot -= f * m_wall[i]->vdot(t); - m_Qdot += f * m_wall[i]->Q(t); + m_vdot -= f * m_wall[i]->expansionRate(); + m_Qdot += f * m_wall[i]->heatRate(); } } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 973cb4eb37..c579e14cc4 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -28,6 +28,7 @@ ReactorNet::~ReactorNet() void ReactorNet::setInitialTime(double time) { m_time = time; + m_initial_time = time; m_integrator_init = false; } diff --git a/src/zeroD/Wall.cpp b/src/zeroD/Wall.cpp index f76fb2f688..8c869ab2a5 100644 --- a/src/zeroD/Wall.cpp +++ b/src/zeroD/Wall.cpp @@ -34,8 +34,17 @@ void WallBase::setArea(double a) { m_surf[1].setArea(a); } +double Wall::velocity() const { + if (m_vf) { + return m_vf->eval(m_time); + } + return 0.; +} + double Wall::vdot(double t) { + warn_deprecated("Wall::vdot", + "To be removed after Cantera 3.0; replaceable by 'expansionRate'."); double rate = m_k * m_area * (m_left->pressure() - m_right->pressure()); if (m_vf) { @@ -44,8 +53,27 @@ double Wall::vdot(double t) return rate; } +double Wall::expansionRate() +{ + double rate = m_k * m_area * (m_left->pressure() - m_right->pressure()); + + if (m_vf) { + rate += m_area * m_vf->eval(m_time); + } + return rate; +} + +double Wall::heatFlux() const { + if (m_qf) { + return m_qf->eval(m_time); + } + return 0.; +} + double Wall::Q(double t) { + warn_deprecated("Wall::Q", + "To be removed after Cantera 3.0; replaceable by 'heatRate'."); double q1 = (m_area * m_rrth) * (m_left->temperature() - m_right->temperature()); if (m_emiss > 0.0) { @@ -60,4 +88,20 @@ double Wall::Q(double t) return q1; } +double Wall::heatRate() +{ + double q1 = (m_area * m_rrth) * + (m_left->temperature() - m_right->temperature()); + if (m_emiss > 0.0) { + double tl = m_left->temperature(); + double tr = m_right->temperature(); + q1 += m_emiss * m_area * StefanBoltz * (tl*tl*tl*tl - tr*tr*tr*tr); + } + + if (m_qf) { + q1 += m_area * m_qf->eval(m_time); + } + return q1; +} + } diff --git a/src/zeroD/flowControllers.cpp b/src/zeroD/flowControllers.cpp index 518a1bf25e..5e74dcd993 100644 --- a/src/zeroD/flowControllers.cpp +++ b/src/zeroD/flowControllers.cpp @@ -31,6 +31,13 @@ void MassFlowController::updateMassFlowRate(double time) m_mdot = std::max(mdot, 0.0); } +void PressureController::setMaster(FlowDevice* master) +{ + warn_deprecated("PressureController::setMaster", + "To be removed after Cantera 3.0; replaced by setPrimary."); + m_primary = master; +} + void PressureController::updateMassFlowRate(double time) { if (!ready()) { @@ -44,8 +51,8 @@ void PressureController::updateMassFlowRate(double time) } else { mdot *= delta_P; } - m_master->updateMassFlowRate(time); - mdot += m_master->massFlowRate(); + m_primary->updateMassFlowRate(time); + mdot += m_primary->massFlowRate(); m_mdot = std::max(mdot, 0.0); } diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 4d3cf0b043..a50428c1a0 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -17,6 +17,7 @@ def make_reactors(self, independent=True, n_reactors=2, T2=300, P2=101325, X2='O2:1.0'): self.net = ct.ReactorNet() + assert self.net.initial_time == 0. self.gas1 = ct.Solution('h2o2.yaml', transport_model=None) self.gas1.TPX = T1, P1, X1 @@ -197,7 +198,8 @@ def test_timestepping(self): self.net.max_time_step = dt_max self.assertEqual(self.net.max_time_step, dt_max) - self.net.set_initial_time(tStart) + self.net.initial_time = tStart + assert self.net.initial_time == tStart self.assertNear(self.net.time, tStart) while t < tEnd: @@ -215,7 +217,7 @@ def test_maxsteps(self): # enough time-steps to reach the endtime max_steps = 10 max_step_size = 1e-07 - self.net.set_initial_time(0) + self.net.initial_time = 0. self.net.max_time_step = max_step_size self.net.max_steps = max_steps with self.assertRaisesRegex( @@ -361,10 +363,10 @@ def test_heat_transfer2(self): w = self.add_wall(U=100, A=0.5) self.assertNear(w.heat_transfer_coeff * w.area * (self.r1.T - self.r2.T), - w.qdot(0)) + w.heat_rate) self.net.advance(1.0) self.assertNear(w.heat_transfer_coeff * w.area * (self.r1.T - self.r2.T), - w.qdot(1.0)) + w.heat_rate) T1b = self.r1.T T2b = self.r2.T @@ -427,19 +429,14 @@ def test_wall_velocity(self): self.add_wall(A=A) - def v(t): - if 0 < t <= 1: - return t - elif 1 <= t <= 2: - return 2 - t - else: - return 0.0 + v = ct.Tabulated1([0.0, 1.0, 2.0], [0.0, 1.0, 0.0]) - self.w.set_velocity(v) + self.w.velocity = v self.net.advance(1.0) - self.assertNear(self.w.vdot(1.0), 1.0 * A, 1e-7) + assert self.w.velocity == approx(v(1.0)) + self.assertNear(self.w.expansion_rate, 1.0 * A, 1e-7) self.net.advance(2.0) - self.assertNear(self.w.vdot(2.0), 0.0, 1e-7) + self.assertNear(self.w.expansion_rate, 0.0, 1e-7) self.assertNear(self.r1.volume, V1 + 1.0 * A, 1e-7) self.assertNear(self.r2.volume, V2 - 1.0 * A, 1e-7) @@ -474,10 +471,12 @@ def test_heat_flux_func(self): V2a = self.r2.volume self.add_wall(A=0.3) - self.w.set_heat_flux(lambda t: 90000 * (1 - t**2) if t <= 1.0 else 0.0) + hfunc = lambda t: 90000 * (1 - t**2) if t <= 1.0 else 0.0 + self.w.heat_flux = hfunc Q = 0.3 * 60000 self.net.advance(1.1) + assert self.w.heat_flux == hfunc(1.1) U1b = self.r1.volume * self.r1.density * self.r1.thermo.u U2b = self.r2.volume * self.r2.density * self.r2.thermo.u @@ -550,7 +549,7 @@ def test_mass_flow_controller_errors(self): self.net.step() with pytest.raises(NotImplementedError): - mfc.set_pressure_function(lambda p: p**2) + mfc.pressure_function = lambda p: p**2 def test_valve1(self): self.make_reactors(P1=10*ct.one_atm, X1='AR:1.0', X2='O2:1.0') @@ -627,7 +626,7 @@ def test_valve3(self): self.net.atol = 1e-20 valve = ct.Valve(self.r1, self.r2) mdot = lambda dP: 5e-3 * np.sqrt(dP) if dP > 0 else 0.0 - valve.set_pressure_function(mdot) + valve.pressure_function = mdot self.assertEqual(valve.valve_coeff, 1.) Y1 = self.r1.Y @@ -656,17 +655,26 @@ def test_valve_timing(self): valve = ct.Valve(self.r1, self.r2) k = 2e-5 valve.valve_coeff = k - valve.set_time_function(lambda t: t>.01) + valve.time_function = lambda t: t > .01 + delta_p = lambda: self.r1.thermo.P - self.r2.thermo.P mdot = lambda: valve.valve_coeff * (self.r1.thermo.P - self.r2.thermo.P) self.net.initialize() - self.assertEqual(valve.mass_flow_rate, 0.0) + assert valve.time_function == 0.0 + assert valve.pressure_function == approx(delta_p()) + assert valve.mass_flow_rate == 0.0 self.net.advance(0.01) - self.assertEqual(valve.mass_flow_rate, 0.0) + assert valve.time_function == 0.0 + assert valve.pressure_function == approx(delta_p()) + assert valve.mass_flow_rate == 0.0 self.net.advance(0.01 + 1e-9) - self.assertNear(valve.mass_flow_rate, mdot()) + assert valve.time_function == 1.0 + assert valve.pressure_function == approx(delta_p()) + assert valve.mass_flow_rate == approx(mdot()) self.net.advance(0.02) - self.assertNear(valve.mass_flow_rate, mdot()) + assert valve.time_function == 1.0 + assert valve.pressure_function == approx(delta_p()) + assert valve.mass_flow_rate == approx(mdot()) def test_valve_errors(self): self.make_reactors() @@ -692,10 +700,10 @@ def test_pressure_controller1(self): mfc = ct.MassFlowController(inlet_reservoir, self.r1) mdot = lambda t: np.exp(-100*(t-0.5)**2) mfc.mass_flow_coeff = 1. - mfc.set_time_function(mdot) + mfc.time_function = mdot pc = ct.PressureController(self.r1, outlet_reservoir) - pc.set_master(mfc) + pc.primary = mfc pc.pressure_coeff = 1e-5 self.assertEqual(pc.pressure_coeff, 1e-5) @@ -717,12 +725,12 @@ def test_pressure_controller2(self): mfc = ct.MassFlowController(inlet_reservoir, self.r1) mdot = lambda t: np.exp(-100*(t-0.5)**2) mfc.mass_flow_coeff = 1. - mfc.set_time_function(mdot) + mfc.time_function = mdot pc = ct.PressureController(self.r1, outlet_reservoir) - pc.set_master(mfc) + pc.primary = mfc pfunc = lambda dp: 1.e-5 * abs(dp)**.5 - pc.set_pressure_function(pfunc) + pc.pressure_function = pfunc self.assertEqual(pc.pressure_coeff, 1.) t = 0 @@ -737,7 +745,7 @@ def test_pressure_controller_errors(self): res = ct.Reservoir(self.gas1) mfc = ct.MassFlowController(res, self.r1, mdot=0.6) - p = ct.PressureController(self.r1, self.r2, master=mfc, K=0.5) + p = ct.PressureController(self.r1, self.r2, primary=mfc, K=0.5) with self.assertRaisesRegex(ct.CanteraError, 'is not ready'): p = ct.PressureController(self.r1, self.r2, K=0.5) @@ -749,14 +757,14 @@ def test_pressure_controller_errors(self): with pytest.raises(NotImplementedError): p = ct.PressureController(self.r1, self.r2) - p.set_time_function(lambda t: t>1.) + p.time_function = lambda t: t>1. def test_set_initial_time(self): self.make_reactors(P1=10*ct.one_atm, X1='AR:1.0', X2='O2:1.0') self.net.rtol = 1e-12 valve = ct.Valve(self.r1, self.r2) - mdot = lambda dP: 5e-3 * np.sqrt(dP) if dP > 0 else 0.0 - valve.set_pressure_function(mdot) + pfunc_a = lambda dP: 5e-3 * np.sqrt(dP) if dP > 0 else 0.0 + valve.pressure_function = pfunc_a t0 = 0.0 tf = t0 + 0.5 @@ -764,20 +772,22 @@ def test_set_initial_time(self): self.assertNear(self.net.time, tf) p1a = self.r1.thermo.P p2a = self.r2.thermo.P + assert valve.pressure_function == approx(pfunc_a(p1a - p2a)) self.make_reactors(P1=10*ct.one_atm, X1='AR:1.0', X2='O2:1.0') self.net.rtol = 1e-12 valve = ct.Valve(self.r1, self.r2) - mdot = lambda dP: 5e-3 * np.sqrt(dP) if dP > 0 else 0.0 - valve.set_pressure_function(mdot) + pfunc_b = lambda dP: 5e-3 * np.sqrt(dP) if dP > 0 else 0.0 + valve.pressure_function = pfunc_b t0 = 0.2 - self.net.set_initial_time(t0) + self.net.initial_time = t0 tf = t0 + 0.5 self.net.advance(tf) self.assertNear(self.net.time, tf) p1b = self.r1.thermo.P p2b = self.r2.thermo.P + assert valve.pressure_function == approx(pfunc_b(p1b - p2b)) self.assertNear(p1a, p1b) self.assertNear(p2a, p2b) @@ -1432,7 +1442,7 @@ def test_reacting(self): while net.distance < 1.0: net.step() i += 1 - assert r.speed * r.density * r.area == pytest.approx(10) + assert r.speed * r.density * r.area == approx(10) stats = net.solver_stats assert stats['steps'] == i @@ -1702,7 +1712,7 @@ def test_reinitialization(self): r.mass_flow_rate = 0.01 r.syncState() - sim.set_initial_time(0.0) + sim.initial_time = 0. sim.advance(0.6) Ygas2 = r.thermo.Y cov2 = rsurf.kinetics.coverages @@ -2473,7 +2483,7 @@ def after_get_state(self, y): def after_update_state(self, y): self.v_wall = y[self.i_wall] - self.walls[0].set_velocity(self.v_wall) + self.walls[0].velocity = self.v_wall def after_eval(self, t, LHS, RHS): # Extra equation is d(v_wall)/dt = k * delta P @@ -2644,12 +2654,13 @@ def after_eval(self,t,LHS,RHS): self.assertNear(add_heat, r1_heat, atol=1e-5) def test_heat_addition(self): - # Applying heat via 'qdot' property should be equivalent to adding it via a wall + # Applying heat via 'heat_rate' property should be equivalent to adding it via + # a wall Qext = 100 Qwall = -66 class HeatedReactor(ct.ExtensibleReactor): def after_eval_walls(self, y): - self.qdot += Qext + self.heat_rate += Qext self.gas.TPX = 300, ct.one_atm, "N2:1.0" r1 = HeatedReactor(self.gas) @@ -2661,7 +2672,7 @@ def after_eval_walls(self, y): net.advance(t) U = r1.thermo.int_energy_mass * r1.mass self.assertNear(U - U0, (Qext + Qwall) * t) - self.assertNear(r1.qdot, Qext + Qwall) + self.assertNear(r1.heat_rate, Qext + Qwall) def test_with_surface(self): phase_defs = """ @@ -2764,8 +2775,8 @@ def after_update_connected(self, do_pressure): def replace_eval_walls(self, t): if self.neighbor: - self.vdot = np.clip(self.p_coeff * (self.P - self.neighbor.P), - -1.7, 1.7) + self.expansion_rate = np.clip( + self.p_coeff * (self.P - self.neighbor.P), -1.7, 1.7) def after_eval(self, t, LHS, RHS): if self.neighbor: @@ -2794,14 +2805,14 @@ def deltaC(): V0 = r1.volume + r2.volume for t in np.linspace(0.01, 0.2, 50): net.advance(t) - self.assertNear(r1.vdot, -r2.vdot) + self.assertNear(r1.expansion_rate, -r2.expansion_rate) self.assertNear(V0, r1.volume + r2.volume) deltaCnow = deltaC() self.assertLess(deltaCnow, deltaCprev) # difference is always decreasing deltaCprev = deltaCnow self.assertArrayNear(M0, r1.mass * r1.thermo.Y + r2.mass * r2.thermo.Y, rtol=2e-8) - states1.append(r1.thermo.state, t=net.time, mass=r1.mass, vdot=r1.vdot) - states2.append(r2.thermo.state, t=net.time, mass=r2.mass, vdot=r2.vdot) + states1.append(r1.thermo.state, t=net.time, mass=r1.mass, vdot=r1.expansion_rate) + states2.append(r2.thermo.state, t=net.time, mass=r2.mass, vdot=r2.expansion_rate) # Regression test values self.assertNear(r1.thermo.P, 151561.15, rtol=1e-6)