diff --git a/doc/sphinx/python/zerodim.rst b/doc/sphinx/python/zerodim.rst index 1737f6e139..c1f02ab346 100644 --- a/doc/sphinx/python/zerodim.rst +++ b/doc/sphinx/python/zerodim.rst @@ -123,6 +123,14 @@ ReactorSurface ^^^^^^^^^^^^^^ .. autoclass:: ReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None) +ExtensibleReactorSurface +^^^^^^^^^^^^^^^^^^^^^^^^ +.. autoclass:: ExtensibleReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None) + +ExtensibleMoleReactorSurface +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. autoclass:: ExtensibleMoleReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None) + Flow Controllers ---------------- diff --git a/include/cantera/base/ct_defs.h b/include/cantera/base/ct_defs.h index bf3e0e3274..2d3052873c 100644 --- a/include/cantera/base/ct_defs.h +++ b/include/cantera/base/ct_defs.h @@ -27,6 +27,7 @@ #include #include #include +#include /** * Namespace for the Cantera kernel. @@ -45,6 +46,7 @@ using std::map; using std::set; using std::function; using std::pair; +using std::span; /** * @defgroup physConstants Physical Constants diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index a546e8f61a..2978e29760 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -274,6 +274,8 @@ class Kinetics return m_start[n] + k; } + size_t speciesOffset(const ThermoPhase& phase) const; + //! Return the name of the kth species in the kinetics manager. /*! * k is an integer from 0 to ktot - 1, where ktot is the number of diff --git a/include/cantera/numerics/EigenSparseJacobian.h b/include/cantera/numerics/EigenSparseJacobian.h index de4232b4f5..39e29e9f08 100644 --- a/include/cantera/numerics/EigenSparseJacobian.h +++ b/include/cantera/numerics/EigenSparseJacobian.h @@ -23,6 +23,11 @@ class EigenSparseJacobian : public SystemJacobian void updatePreconditioner() override; void updateTransient(double rdt, int* mask) override; + //! Set all Jacobian elements. Replaces the existing elements. + void setFromTriplets(const vector>& trips) { + m_jac_trips = trips; + } + //! Return underlying Jacobian matrix //! @ingroup derivGroup Eigen::SparseMatrix jacobian(); diff --git a/include/cantera/zeroD/ConstPressureMoleReactor.h b/include/cantera/zeroD/ConstPressureMoleReactor.h index 5031e48899..0ea9c0ce20 100644 --- a/include/cantera/zeroD/ConstPressureMoleReactor.h +++ b/include/cantera/zeroD/ConstPressureMoleReactor.h @@ -20,16 +20,15 @@ namespace Cantera class ConstPressureMoleReactor : public MoleReactor { public: - using MoleReactor::MoleReactor; // inherit constructors + ConstPressureMoleReactor(shared_ptr sol, const string& name="(none)"); + ConstPressureMoleReactor(shared_ptr sol, bool clone, + const string& name="(none)"); string type() const override { return "ConstPressureMoleReactor"; }; void getState(double* y) override; - - void initialize(double t0=0.0) override; - void eval(double t, double* LHS, double* RHS) override; vector steadyConstraints() const override { diff --git a/include/cantera/zeroD/ConstPressureReactor.h b/include/cantera/zeroD/ConstPressureReactor.h index a443cc6a45..9566289b74 100644 --- a/include/cantera/zeroD/ConstPressureReactor.h +++ b/include/cantera/zeroD/ConstPressureReactor.h @@ -23,15 +23,15 @@ namespace Cantera class ConstPressureReactor : public Reactor { public: - using Reactor::Reactor; // inherit constructors + ConstPressureReactor(shared_ptr sol, const string& name="(none)"); + ConstPressureReactor(shared_ptr sol, bool clone, + const string& name="(none)"); string type() const override { return "ConstPressureReactor"; } void getState(double* y) override; - - void initialize(double t0=0.0) override; void eval(double t, double* LHS, double* RHS) override; vector steadyConstraints() const override; diff --git a/include/cantera/zeroD/FlowReactor.h b/include/cantera/zeroD/FlowReactor.h index dfb97e5629..ea2f4a481a 100644 --- a/include/cantera/zeroD/FlowReactor.h +++ b/include/cantera/zeroD/FlowReactor.h @@ -16,7 +16,8 @@ namespace Cantera class FlowReactor : public IdealGasReactor { public: - using IdealGasReactor::IdealGasReactor; // inherit constructors + FlowReactor(shared_ptr sol, const string& name="(none)"); + FlowReactor(shared_ptr sol, bool clone, const string& name="(none)"); string type() const override { return "FlowReactor"; @@ -36,8 +37,6 @@ class FlowReactor : public IdealGasReactor } void getStateDae(double* y, double* ydot) override; - void initialize(double t0=0.0) override; - void syncState() override; void updateState(double* y) override; //! Not implemented; FlowReactor implements evalDae() instead. @@ -74,102 +73,26 @@ class FlowReactor : public IdealGasReactor //! Sets the area of the reactor [m²] void setArea(double area) override; - //! The ratio of the reactor's surface area to volume ratio [m^-1] - //! @note If the surface area to volume ratio is unspecified by the user, - //! this will be calculated assuming the reactor is a cylinder. - double surfaceAreaToVolumeRatio() const; - - //! Set the reactor's surface area to volume ratio [m^-1] - void setSurfaceAreaToVolumeRatio(double sa_to_vol) { - m_sa_to_vol = sa_to_vol; - } - - //! Get the steady state tolerances used to determine the initial state for - //! surface coverages - double inletSurfaceAtol() const { - return m_ss_atol; - } - - //! Set the steady state tolerances used to determine the initial state for - //! surface coverages - void setInletSurfaceAtol(double atol) { - m_ss_atol = atol; - } - - //! Get the steady state tolerances used to determine the initial state for - //! surface coverages - double inletSurfaceRtol() const { - return m_ss_rtol; - } - - //! Set the steady state tolerances used to determine the initial state for - //! surface coverages - void setInletSurfaceRtol(double rtol) { - m_ss_rtol = rtol; - } - - //! Get the steady state tolerances used to determine the initial state for - //! surface coverages - double inletSurfaceMaxSteps() const { - return m_max_ss_steps; - } - - //! Set the steady state tolerances used to determine the initial state for - //! surface coverages - void setInletSurfaceMaxSteps(int max_steps) { - m_max_ss_steps = max_steps; - } - - //! Get the steady state tolerances used to determine the initial state for - //! surface coverages - double inletSurfaceMaxErrorFailures() const { - return m_max_ss_error_fails; - } - - //! Set the steady state tolerances used to determine the initial state for - //! surface coverages - void setInletSurfaceMaxErrorFailures(int max_fails) { - m_max_ss_error_fails = max_fails; - } - //! Return the index in the solution vector for this reactor of the component named //! *nm*. Possible values for *nm* are "density", "speed", "pressure", - //! "temperature", the name of a homogeneous phase species, or the name of a surface - //! species. + //! "temperature" or the name of a homogeneous phase species. size_t componentIndex(const string& nm) const override; string componentName(size_t k) override; - void updateSurfaceState(double* y) override; - protected: //! Density [kg/m^3]. First component of the state vector. double m_rho = NAN; //! Axial velocity [m/s]. Second component of the state vector. double m_u = -1.0; - //! Pressure [Pa]. Third component of the state vector. - double m_P = NAN; - //! Temperature [K]. Fourth component of the state vector. - double m_T = NAN; //! offset to the species equations const size_t m_offset_Y = 4; //! reactor area [m^2] double m_area = 1.0; - //! reactor surface area to volume ratio [m^-1] - double m_sa_to_vol = -1.0; - //! temporary storage for surface species production rates - vector m_sdot_temp; //! temporary storage for species partial molar enthalpies vector m_hk; - //! steady-state relative tolerance, used to determine initial surface coverages - double m_ss_rtol = 1e-7; - //! steady-state absolute tolerance, used to determine initial surface coverages - double m_ss_atol = 1e-14; - //! maximum number of steady-state coverage integrator-steps - int m_max_ss_steps = 20000; - //! maximum number of steady-state integrator error test failures - int m_max_ss_error_fails = 10; }; + } #endif diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 27de6d57c4..1fbb1b1c58 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -33,13 +33,14 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor void eval(double t, double* LHS, double* RHS) override; void updateState(double* y) override; + void getJacobianScalingFactors(double& f_species, double* f_energy) override; //! Calculate an approximate Jacobian to accelerate preconditioned solvers //! Neglects derivatives with respect to mole fractions that would generate a //! fully-dense Jacobian. Currently also neglects terms related to interactions //! between reactors, for example via inlets and outlets. - Eigen::SparseMatrix jacobian() override; + void getJacobianElements(vector>& trips) override; bool preconditionerSupported() const override { return true; }; @@ -50,6 +51,7 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor protected: vector m_hk; //!< Species molar enthalpies + double m_TotalCp; //!< Total heat capacity (@f$ m c_p @f$) [J/K] }; } diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index 9a2fceece4..3347eaa6f8 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -45,12 +45,14 @@ class IdealGasMoleReactor : public MoleReactor //! Neglects derivatives with respect to mole fractions that would generate a //! fully-dense Jacobian. Currently, also neglects terms related to interactions //! between reactors, for example via inlets and outlets. - Eigen::SparseMatrix jacobian() override; + void getJacobianElements(vector>& trips) override; + void getJacobianScalingFactors(double& f_species, double* f_energy) override; bool preconditionerSupported() const override {return true;}; protected: vector m_uk; //!< Species molar internal energies + double m_TotalCv; //!< Total heat capacity (@f$ m c_v @f$) [J/K] }; } diff --git a/include/cantera/zeroD/MoleReactor.h b/include/cantera/zeroD/MoleReactor.h index e15e07695d..4a9719751b 100644 --- a/include/cantera/zeroD/MoleReactor.h +++ b/include/cantera/zeroD/MoleReactor.h @@ -20,20 +20,16 @@ namespace Cantera class MoleReactor : public Reactor { public: - using Reactor::Reactor; // inherit constructors + MoleReactor(shared_ptr sol, const string& name="(none)"); + MoleReactor(shared_ptr sol, bool clone, const string& name="(none)"); string type() const override { return "MoleReactor"; } - void initialize(double t0=0.0) override; - void getState(double* y) override; - void updateState(double* y) override; - void eval(double t, double* LHS, double* RHS) override; - size_t componentIndex(const string& nm) const override; string componentName(size_t k) override; double upperBound(size_t k) const override; @@ -41,11 +37,6 @@ class MoleReactor : public Reactor void resetBadValues(double* y) override; protected: - //! For each surface in the reactor, update vector of triplets with all relevant - //! surface jacobian derivatives of species with respect to species - //! which are appropriately offset to align with the reactor's state vector. - virtual void addSurfaceJacobian(vector> &triplets); - //! Get moles of the system from mass fractions stored by thermo object //! @param y vector for moles to be put into void getMoles(double* y); @@ -54,12 +45,6 @@ class MoleReactor : public Reactor //! @param y vector of moles of the system void setMassFromMoles(double* y); - void evalSurfaces(double* LHS, double* RHS, double* sdot) override; - - void updateSurfaceState(double* y) override; - - void getSurfaceInitialConditions(double* y) override; - //! const value for the species start index const size_t m_sidx = 2; }; diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 8dcaf5853d..11004fac09 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -60,12 +60,6 @@ class Reactor : public ReactorBase return true; } - //! Indicates whether the governing equations for this reactor are functions of time - //! or a spatial variable. All reactors in a network must have the same value. - virtual bool timeIsIndependent() const { - return true; - } - void setInitialVolume(double vol) override { m_vol = vol; } @@ -86,98 +80,23 @@ class Reactor : public ReactorBase return m_energy; } - //! Number of equations (state variables) for this reactor - size_t neq() { - if (!m_nv) { - initialize(); - } - return m_nv; - } - - //! Get the the current state of the reactor. - /*! - * @param[out] y state vector representing the initial state of the reactor - */ - virtual void getState(double* y); - - //! Get the current state and derivative vector of the reactor for a DAE solver - /*! - * @param[out] y state vector representing the initial state of the reactor - * @param[out] ydot state vector representing the initial derivatives of the - * reactor - */ - virtual void getStateDae(double* y, double* ydot) { - throw NotImplementedError("Reactor::getStateDae(y, ydot)"); - } - + void getState(double* y) override; void initialize(double t0=0.0) override; - - //! Evaluate the reactor governing equations. Called by ReactorNet::eval. - //! @param[in] t time. - //! @param[out] LHS pointer to start of vector of left-hand side - //! coefficients for governing equations, length m_nv, default values 1 - //! @param[out] RHS pointer to start of vector of right-hand side - //! coefficients for governing equations, length m_nv, default values 0 - virtual void eval(double t, double* LHS, double* RHS); - - /** - * Evaluate the reactor governing equations. Called by ReactorNet::eval. - * @param[in] t time. - * @param[in] y solution vector, length neq() - * @param[in] ydot rate of change of solution vector, length neq() - * @param[out] residual residuals vector, length neq() - */ - virtual void evalDae(double t, double* y, double* ydot, double* residual) { - throw NotImplementedError("Reactor::evalDae"); - } - - //! Given a vector of length neq(), mark which variables should be - //! considered algebraic constraints - virtual void getConstraints(double* constraints) { - throw NotImplementedError("Reactor::getConstraints"); - } - - //! Get the indices of equations that are algebraic constraints when solving the - //! steady-state problem. - //! - //! @warning This method is an experimental part of the %Cantera API and may be - //! changed or removed without notice. - //! @since New in %Cantera 3.2. - virtual vector steadyConstraints() const; - - //! Set the state of the reactor to correspond to the state vector *y*. - virtual void updateState(double* y); - - //! Number of sensitivity parameters associated with this reactor. - //! (including walls) - size_t nSensParams() const override; - + void eval(double t, double* LHS, double* RHS) override; + vector steadyConstraints() const override; + void updateState(double* y) override; void addSensitivityReaction(size_t rxn) override; - - //! Add a sensitivity parameter associated with the enthalpy formation of - //! species *k* (in the homogeneous phase) - virtual void addSensitivitySpeciesEnthalpy(size_t k); + void addSensitivitySpeciesEnthalpy(size_t k) override; //! Return the index in the solution vector for this reactor of the //! component named *nm*. Possible values for *nm* are "mass", "volume", //! "int_energy", the name of a homogeneous phase species, or the name of a //! surface species. - virtual size_t componentIndex(const string& nm) const; - - //! Return the name of the solution component with index *i*. - //! @see componentIndex() - virtual string componentName(size_t k); - - //! Get the upper bound on the k-th component of the local state vector. - virtual double upperBound(size_t k) const; - - //! Get the lower bound on the k-th component of the local state vector. - virtual double lowerBound(size_t k) const; - - //! Reset physically or mathematically problematic values, such as negative species - //! concentrations. - //! @param[inout] y current state vector, to be updated; length neq() - virtual void resetBadValues(double* y); + size_t componentIndex(const string& nm) const override; + string componentName(size_t k) override; + double upperBound(size_t k) const override; + double lowerBound(size_t k) const override; + void resetBadValues(double* y) override; //! Set absolute step size limits during advance //! @param limits array of step size limits with length neq @@ -199,18 +118,6 @@ class Reactor : public ReactorBase //! @param limit value for step size limit void setAdvanceLimit(const string& nm, const double limit); - //! Calculate the Jacobian of a specific Reactor specialization. - //! @warning Depending on the particular implementation, this may return an - //! approximate Jacobian intended only for use in forming a preconditioner for - //! iterative solvers. - //! @ingroup derivGroup - //! - //! @warning This method is an experimental part of the %Cantera - //! API and may be changed or removed without notice. - virtual Eigen::SparseMatrix jacobian() { - throw NotImplementedError("Reactor::jacobian"); - } - //! Calculate the reactor-specific Jacobian using a finite difference method. //! //! This method is used only for informational purposes. Jacobian calculations @@ -220,84 +127,36 @@ class Reactor : public ReactorBase //! API and may be changed or removed without notice. Eigen::SparseMatrix finiteDifferenceJacobian(); - //! Use this to set the kinetics objects derivative settings - virtual void setDerivativeSettings(AnyMap& settings); - - //! Set reaction rate multipliers based on the sensitivity variables in - //! *params*. - virtual void applySensitivity(double* params); - - //! Reset the reaction rate multipliers - virtual void resetSensitivity(double* params); + void setDerivativeSettings(AnyMap& settings) override; - //! Return a false if preconditioning is not supported or true otherwise. - //! - //! @warning This method is an experimental part of the %Cantera - //! API and may be changed or removed without notice. - //! - //! @since New in %Cantera 3.0 - //! - virtual bool preconditionerSupported() const {return false;}; + void applySensitivity(double* params) override; + void resetSensitivity(double* params) override; protected: - //! Return the index in the solution vector for this reactor of the species - //! named *nm*, in either the homogeneous phase or a surface phase, relative - //! to the start of the species terms. Used to implement componentIndex for - //! specific reactor implementations. - virtual size_t speciesIndex(const string& nm) const; + //! Update #m_sdot to reflect current production rates of bulk phase species due to + //! reactions on adjacent surfaces. + void updateSurfaceProductionRates(); //! Evaluate terms related to Walls. Calculates #m_vdot and #m_Qdot based on //! wall movement and heat transfer. //! @param t the current time - virtual void evalWalls(double t); - - //! Evaluate terms related to surface reactions. - //! @param[out] LHS Multiplicative factor on the left hand side of ODE for surface - //! species coverages - //! @param[out] RHS Right hand side of ODE for surface species coverages - //! @param[out] sdot array of production rates of bulk phase species on surfaces - //! [kmol/s] - virtual void evalSurfaces(double* LHS, double* RHS, double* sdot); - - virtual void evalSurfaces(double* RHS, double* sdot); - - //! Update the state of SurfPhase objects attached to this reactor - virtual void updateSurfaceState(double* y); - - //! 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 - //! or an independent variable. - virtual void updateConnected(bool updatePressure); - - //! Get initial conditions for SurfPhase objects attached to this reactor - virtual void getSurfaceInitialConditions(double* y); + void evalWalls(double t) override; //! Pointer to the homogeneous Kinetics object that handles the reactions Kinetics* m_kin = nullptr; double m_vdot = 0.0; //!< net rate of volume change from moving walls [m^3/s] - double m_Qdot = 0.0; //!< net heat transfer into the reactor, through walls [W] + vector m_wdot; //!< Species net molar production rates + vector m_uk; //!< Species molar internal energies - vector m_work; - - //! Production rates of gas phase species on surfaces [kmol/s] + //! Total production rate of bulk phase species on surfaces [kmol/s] vector m_sdot; - vector m_wdot; //!< Species net molar production rates - vector m_uk; //!< Species molar internal energies bool m_chem = false; bool m_energy = true; - size_t m_nv = 0; - size_t m_nv_surf; //!!< Number of variables associated with reactor surfaces vector m_advancelimits; //!< Advance step limit - - //! Vector of triplets representing the jacobian - vector> m_jac_trips; }; } diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 0fed7946eb..fb434b416a 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -8,6 +8,7 @@ #include "cantera/base/global.h" #include "cantera/base/ctexceptions.h" +#include "cantera/numerics/eigen_sparse.h" namespace Cantera { @@ -25,6 +26,7 @@ class ReactorSurface; class Kinetics; class ThermoPhase; class Solution; +class AnyMap; enum class SensParameterType { reaction, @@ -45,7 +47,7 @@ struct SensitivityParameter * defined. * @ingroup reactorGroup */ -class ReactorBase +class ReactorBase : public std::enable_shared_from_this { public: //! Instantiate a ReactorBase object with Solution contents. @@ -96,6 +98,17 @@ class ReactorBase //! @since New in %Cantera 3.2 shared_ptr phase() const { return m_solution; } + //! Indicates whether the governing equations for this reactor are functions of time + //! or a spatial variable. All reactors in a network must have the same value. + virtual bool timeIsIndependent() const { + return true; + } + + //! Number of equations (state variables) for this reactor + size_t neq() { + return m_nv; + } + //! @name Methods to set up a simulation //! @{ @@ -112,6 +125,12 @@ class ReactorBase "Area is undefined for reactors of type '{}'.", type()); } + //! Evaluate contributions from walls connected to this reactor + virtual void evalWalls(double t) { + throw NotImplementedError("ReactorBase::evalWalls", + "Wall evaluation is undefined for reactors of type '{}'.", type()); + } + //! Set an area associated with a reactor [m²]. //! Examples: surface area of ReactorSurface or cross section area of FlowReactor. virtual void setArea(double a) { @@ -200,6 +219,18 @@ class ReactorBase return m_surfaces.size(); } + //! Production rates on surfaces. + //! + //! For bulk reactors, this contains the total production rates [kmol/s] of bulk + //! phase species due to reactions on all adjacent species. + //! + //! For surfaces, this contains the production rates [kmol/m²/s] of species on the + //! surface and all adjacent phases, in the order defined by the InterfaceKinetics + //! object. + span surfaceProductionRates() const { + return m_sdot; + } + /** * Initialize the reactor. Called automatically by ReactorNet::initialize. */ @@ -207,17 +238,166 @@ class ReactorBase throw NotImplementedError("ReactorBase::initialize"); } - //! @} + //! Get the current state of the reactor. + /*! + * @param[out] y state vector representing the initial state of the reactor + */ + virtual void getState(double* y) { + throw NotImplementedError("ReactorBase::getState"); + } + + //! Get the current state and derivative vector of the reactor for a DAE solver + /*! + * @param[out] y state vector representing the initial state of the reactor + * @param[out] ydot state vector representing the initial derivatives of the + * reactor + */ + virtual void getStateDae(double* y, double* ydot) { + throw NotImplementedError("ReactorBase::getStateDae(y, ydot)"); + } + + //! Evaluate the reactor governing equations. Called by ReactorNet::eval. + //! @param[in] t time. + //! @param[out] LHS pointer to start of vector of left-hand side + //! coefficients for governing equations, length m_nv, default values 1 + //! @param[out] RHS pointer to start of vector of right-hand side + //! coefficients for governing equations, length m_nv, default values 0 + virtual void eval(double t, double* LHS, double* RHS) { + throw NotImplementedError("ReactorBase::eval"); + } + + /** + * Evaluate the reactor governing equations. Called by ReactorNet::eval. + * @param[in] t time. + * @param[in] y solution vector, length neq() + * @param[in] ydot rate of change of solution vector, length neq() + * @param[out] residual residuals vector, length neq() + */ + virtual void evalDae(double t, double* y, double* ydot, double* residual) { + throw NotImplementedError("ReactorBase::evalDae"); + } + + //! Given a vector of length neq(), mark which variables should be + //! considered algebraic constraints + virtual void getConstraints(double* constraints) { + throw NotImplementedError("ReactorBase::getConstraints"); + } + + //! Get the indices of equations that are algebraic constraints when solving the + //! steady-state problem. + //! + //! @warning This method is an experimental part of the %Cantera API and may be + //! changed or removed without notice. + //! @since New in %Cantera 3.2. + virtual vector steadyConstraints() const { + throw NotImplementedError("ReactorBase::steadyConstraints"); + } + + //! Set the state of the reactor to correspond to the state vector *y*. + virtual void updateState(double* y) { + throw NotImplementedError("ReactorBase::updateState"); + } + + //! Add a sensitivity parameter associated with the enthalpy formation of + //! species *k*. + virtual void addSensitivitySpeciesEnthalpy(size_t k) { + throw NotImplementedError("ReactorBase::addSensitivitySpeciesEnthalpy"); + } + + //! Return the index in the solution vector for this reactor of the + //! component named *nm*. + virtual size_t componentIndex(const string& nm) const { + throw NotImplementedError("ReactorBase::componentIndex"); + } + + //! Return the name of the solution component with index *i*. + //! @see componentIndex() + virtual string componentName(size_t k) { + throw NotImplementedError("ReactorBase::componentName"); + } + + //! Get the upper bound on the k-th component of the local state vector. + virtual double upperBound(size_t k) const { + throw NotImplementedError("ReactorBase::upperBound"); + } + + //! Get the lower bound on the k-th component of the local state vector. + virtual double lowerBound(size_t k) const { + throw NotImplementedError("ReactorBase::lowerBound"); + } - //! Set the state of the Phase object associated with this reactor to the - //! reactor's current state. - virtual void restoreState(); + //! Reset physically or mathematically problematic values, such as negative species + //! concentrations. + //! @param[inout] y current state vector, to be updated; length neq() + virtual void resetBadValues(double* y) { + throw NotImplementedError("ReactorBase::resetBadValues"); + } + + //! Get Jacobian elements for this reactor within the full reactor network. + //! + //! Indices within `trips` are global indices within the full reactor network. The + //! reactor is responsible for providing all elements of the Jacobian in the rows + //! corresponding to its state variables, that is, all derivatives of its state + //! variables with respect to all state variables in the network. + //! + //! @warning This method is an experimental part of the %Cantera API and may be + //! changed or removed without notice. + virtual void getJacobianElements(vector>& trips) {}; + + //! Calculate the Jacobian of a specific reactor specialization. + //! @warning Depending on the particular implementation, this may return an + //! approximate Jacobian intended only for use in forming a preconditioner for + //! iterative solvers. + //! @ingroup derivGroup + //! + //! @warning This method is an experimental part of the %Cantera + //! API and may be changed or removed without notice. + virtual Eigen::SparseMatrix jacobian(); + + //! Use this to set the kinetics objects derivative settings + virtual void setDerivativeSettings(AnyMap& settings) { + throw NotImplementedError("ReactorBase::setDerivativeSettings"); + } + + //! Set reaction rate multipliers based on the sensitivity variables in + //! *params*. + virtual void applySensitivity(double* params) { + throw NotImplementedError("ReactorBase::applySensitivity"); + } + + //! Reset the reaction rate multipliers + virtual void resetSensitivity(double* params) { + throw NotImplementedError("ReactorBase::resetSensitivity"); + } + + //! Return a false if preconditioning is not supported or true otherwise. + //! + //! @warning This method is an experimental part of the %Cantera + //! API and may be changed or removed without notice. + //! + //! @since New in %Cantera 3.0 + //! + virtual bool preconditionerSupported() const { return false; }; + + //! @} //! Set the state of the reactor to the associated ThermoPhase object. - //! This method is the inverse of restoreState() and will trigger integrator - //! reinitialization. + //! This method will trigger integrator reinitialization. + //! @deprecated To be removed after %Cantera 4.0. Use ReactorNet::reinitialize to + //! indicate a change in state that requires integrator reinitialization. virtual void syncState(); + //! Update state information needed by connected reactors, flow devices, and walls. + //! + //! Called from updateState() for normal reactor types, and from + //! ReactorNet::updateState for Reservoir. + //! + //! @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); + //! Return the residence time (s) of the contents of this reactor, based //! on the outlet mass flow rates and the mass of the reactor contents. double residenceTime(); @@ -234,22 +414,10 @@ class ReactorBase } //! Returns the current density (kg/m^3) of the reactor's contents. - double density() const { - if (m_state.empty()) { - throw CanteraError("ReactorBase::density", - "Reactor state empty and/or contents not defined."); - } - return m_state[1]; - } + double density() const; //! Returns the current temperature (K) of the reactor's contents. - double temperature() const { - if (m_state.empty()) { - throw CanteraError("ReactorBase::temperature", - "Reactor state empty and/or contents not defined."); - } - return m_state[0]; - } + double temperature() const; //! Returns the current enthalpy (J/kg) of the reactor's contents. double enthalpy_mass() const { @@ -267,22 +435,10 @@ class ReactorBase } //! Return the vector of species mass fractions. - const double* massFractions() const { - if (m_state.empty()) { - throw CanteraError("ReactorBase::massFractions", - "Reactor state empty and/or contents not defined."); - } - return m_state.data() + 2; - } + const double* massFractions() const; //! Return the mass fraction of the *k*-th species. - double massFraction(size_t k) const { - if (m_state.empty()) { - throw CanteraError("ReactorBase::massFraction", - "Reactor state empty and/or contents not defined."); - } - return m_state[k+2]; - } + double massFraction(size_t k) const; //! @} @@ -292,6 +448,32 @@ class ReactorBase //! Set the ReactorNet that this reactor belongs to. void setNetwork(ReactorNet* net); + //! Get the starting offset for this reactor's state variables within the global + //! state vector of the ReactorNet. + size_t offset() const { return m_offset; } + + //! Set the starting offset for this reactor's state variables within the global + //! state vector of the ReactorNet. + void setOffset(size_t offset) { m_offset = offset; } + + //! Offset of the first species in the local state vector + size_t speciesOffset() const { + return m_nv - m_nsp; + } + + //! Get scaling factors for the Jacobian matrix terms proportional to + //! @f$ d\dot{n}_k/dC_j @f$. + //! + //! Used to determine contribution of surface phases to the Jacobian. + //! + //! @param f_species Scaling factor for derivatives appearing in the species + //! equations. Equal to $1/V$. + //! @param f_energy Scaling factor for each species term appearing in the energy + //! equation. + virtual void getJacobianScalingFactors(double& f_species, double* f_energy) { + throw NotImplementedError("ReactorBase::getJacobianScalingFactors"); + } + //! Add a sensitivity parameter associated with the reaction number *rxn* virtual void addSensitivityReaction(size_t rxn) { throw NotImplementedError("ReactorBase::addSensitivityReaction"); @@ -318,15 +500,17 @@ class ReactorBase vector m_wall; vector m_surfaces; + vector m_sdot; //!< species production rates on surfaces //! Vector of length nWalls(), indicating whether this reactor is on the left (0) //! or right (1) of each wall. vector m_lr; string m_name; //!< Reactor name. bool m_defaultNameSet = false; //!< `true` if default name has been previously set. + size_t m_nv = 0; //!< Number of state variables for this reactor - //! The ReactorNet that this reactor is part of - ReactorNet* m_net = nullptr; + ReactorNet* m_net = nullptr; //!< The ReactorNet that this reactor is part of + size_t m_offset = 0; //!< Offset into global ReactorNet state vector //! Composite thermo/kinetics/transport handler shared_ptr m_solution; diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h index f2556a7871..9574ea6c80 100644 --- a/include/cantera/zeroD/ReactorDelegator.h +++ b/include/cantera/zeroD/ReactorDelegator.h @@ -45,12 +45,8 @@ class ReactorAccessor //! delegate for Reactor::evalWalls(). 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; - - //! Set the state of the thermo object for surface *n* to correspond to the - //! state of that surface - virtual void restoreSurfaceState(size_t n) = 0; + //! @copydoc ReactorBase::surfaceProductionRates + virtual span surfaceProductionRates() = 0; }; //! Delegate methods of the Reactor class to external functions @@ -59,22 +55,15 @@ template class ReactorDelegator : public Delegator, public R, public ReactorAccessor { public: - ReactorDelegator(shared_ptr phase, bool clone, const string& name="(none)") - : R(phase, clone, name) + template + ReactorDelegator(Args&&... args) + : R(std::forward(args)...) { install("initialize", m_initialize, [this](double t0) { R::initialize(t0); }); - install("syncState", m_syncState, [this]() { R::syncState(); }); install("getState", m_getState, [this](std::array sizes, double* y) { R::getState(y); }); install("updateState", m_updateState, [this](std::array sizes, double* y) { R::updateState(y); }); - install("updateSurfaceState", m_updateSurfaceState, - [this](std::array sizes, double* y) { R::updateSurfaceState(y); }); - install("getSurfaceInitialConditions", m_getSurfaceInitialConditions, - [this](std::array sizes, double* y) { - R::getSurfaceInitialConditions(y); - } - ); install("updateConnected", m_updateConnected, [this](bool updatePressure) { R::updateConnected(updatePressure); }); install("eval", m_eval, @@ -82,18 +71,13 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor R::eval(t, LHS, RHS); } ); - install("evalWalls", m_evalWalls, [this](double t) { R::evalWalls(t); }); - install("evalSurfaces", m_evalSurfaces, - [this](std::array sizes, double* LHS, double* RHS, double* sdot) { - R::evalSurfaces(LHS, RHS, sdot); - } - ); + if constexpr (std::is_base_of::value) { + install("evalWalls", m_evalWalls, [this](double t) { R::evalWalls(t); }); + } install("componentName", m_componentName, [this](size_t k) { return R::componentName(k); }); install("componentIndex", m_componentIndex, [this](const string& nm) { return R::componentIndex(nm); }); - install("speciesIndex", m_speciesIndex, - [this](const string& nm) { return R::speciesIndex(nm); }); } // Overrides of Reactor methods @@ -106,10 +90,6 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor m_initialize(t0); } - void syncState() override { - m_syncState(); - } - void getState(double* y) override { std::array sizes{R::neq()}; m_getState(sizes, y); @@ -120,16 +100,6 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor m_updateState(sizes, y); } - void updateSurfaceState(double* y) override { - std::array sizes{R::m_nv_surf}; - m_updateSurfaceState(sizes, y); - } - - void getSurfaceInitialConditions(double* y) override { - std::array sizes{R::m_nv_surf}; - m_getSurfaceInitialConditions(sizes, y); - } - void updateConnected(bool updatePressure) override { m_updateConnected(updatePressure); } @@ -140,12 +110,11 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor } void evalWalls(double t) override { - m_evalWalls(t); - } - - void evalSurfaces(double* LHS, double* RHS, double* sdot) override { - std::array sizes{R::m_nv_surf, R::m_nv_surf, R::m_nsp}; - m_evalSurfaces(sizes, LHS, RHS, sdot); + if constexpr (std::is_base_of::value) { + m_evalWalls(t); + } else { + ReactorBase::evalWalls(t); + } } string componentName(size_t k) override { @@ -156,10 +125,6 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor return m_componentIndex(nm); } - size_t speciesIndex(const string& nm) const override { - return m_speciesIndex(nm); - } - // Public access to protected Reactor variables needed by derived classes void setNEq(size_t n) override { @@ -167,43 +132,54 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor } double expansionRate() const override { - return R::m_vdot; + if constexpr (std::is_base_of::value) { + return R::m_vdot; + } else { + throw NotImplementedError("ReactorDelegator::expansionRate", + "Expansion rate is undefined for reactors of type '{}'.", type()); + } } void setExpansionRate(double v) override { - R::m_vdot = v; + if constexpr (std::is_base_of::value) { + R::m_vdot = v; + } else { + throw NotImplementedError("ReactorDelegator::setExpansionRate", + "Expansion rate is undefined for reactors of type '{}'.", type()); + } } double heatRate() const override { - return R::m_Qdot; + if constexpr (std::is_base_of::value) { + return R::m_Qdot; + } else { + throw NotImplementedError("ReactorDelegator::heatRate", + "Heat rate is undefined for reactors of type '{}'.", type()); + } } void setHeatRate(double q) override { - R::m_Qdot = q; - } - - void restoreThermoState() override { - R::m_thermo->restoreState(R::m_state); + if constexpr (std::is_base_of::value) { + R::m_Qdot = q; + } else { + throw NotImplementedError("ReactorDelegator::setHeatRate", + "Heat rate is undefined for reactors of type '{}'.", type()); + } } - void restoreSurfaceState(size_t n) override { - R::m_surfaces.at(n)->syncState(); + span surfaceProductionRates() override { + return R::m_sdot; } private: function m_initialize; - function m_syncState; function, double*)> m_getState; function, double*)> m_updateState; - function, double*)> m_updateSurfaceState; - function, double*)> m_getSurfaceInitialConditions; function m_updateConnected; function, double, double*, double*)> m_eval; function m_evalWalls; - function, double*, double*, double*)> m_evalSurfaces; function m_componentName; function m_componentIndex; - function m_speciesIndex; }; } diff --git a/include/cantera/zeroD/ReactorFactory.h b/include/cantera/zeroD/ReactorFactory.h index 7b1c6bee5a..b64fe935bd 100644 --- a/include/cantera/zeroD/ReactorFactory.h +++ b/include/cantera/zeroD/ReactorFactory.h @@ -35,6 +35,30 @@ class ReactorFactory : public Factory, bool, c ReactorFactory(); }; +//! Factory class to create ReactorSurface objects +//! +//! This class is mainly used via the newReactorSurface() function, for example: +//! +//! ```cpp +//! shared_ptr rsurf1 = newReactorSurface( +//! surf_phase, {reactor1, reactor2}); +//! @tem +//! ``` +class ReactorSurfaceFactory : public Factory /* surface */, + const vector>& /* adjacent reactors */, + bool /* clone */, const string& /* name */> +{ +public: + static ReactorSurfaceFactory* factory(); + void deleteFactory() override; + +private: + static ReactorSurfaceFactory* s_factory; + static std::mutex s_mutex; + ReactorSurfaceFactory(); +}; + //! @defgroup reactorGroup Reactors //! Zero-dimensional objects representing stirred reactors. //! Reactors simulate time-dependent behavior considering gas-phase chemistry. @@ -87,6 +111,10 @@ shared_ptr newReservoir( //! Create a ReactorSurface object with the specified contents and adjacent reactors //! participating in surface reactions. +//! +//! The type of the reactor surface is automatically determined based on the types of +//! the adjacent reactors. +//! //! @param phase Solution (Interface) object to model the thermodynamic properties //! and reactions occurring in the reactor //! @param reactors List of Reactors adjacent to this surface, whose contents @@ -100,6 +128,23 @@ shared_ptr newReactorSurface( shared_ptr phase, const vector>& reactors, bool clone=true, const string& name="(none)"); +//! Create a ReactorSurface object with the specified contents and adjacent reactors +//! participating in surface reactions. +//! @param model Type of ReactorSurface to be created. +//! @param phase Solution (Interface) object to model the thermodynamic properties +//! and reactions occurring in the reactor +//! @param reactors List of Reactors adjacent to this surface, whose contents +//! participate in reactions occurring on this surface. +//! @param clone Determines whether to clone `sol` so that the internal state of this +//! surface is independent of the original Interface object and any Solution objects +//! used by other reactors in the network except those in the `reactors` list. +//! @param name Name of the reactor surface. +//! @since New in %Cantera 4.0. +shared_ptr newReactorSurface( + const string& model, shared_ptr phase, + const vector>& reactors, + bool clone=true, const string& name="(none)"); + //! @} } diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 2a885a3406..b4736e0668 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -187,7 +187,7 @@ class ReactorNet : public FuncEval //! indices are determined by the order in which the reactors were added //! to the reactor network. Reactor& reactor(int n) { - return *m_reactors[n]; + return *m_bulkReactors[n]; } //! Returns `true` if verbose logging output is enabled. @@ -392,7 +392,12 @@ class ReactorNet : public FuncEval //! and deliberately not exposed in external interfaces. virtual int lastOrder() const; - vector m_reactors; + vector> m_reactors; + vector m_bulkReactors; + vector m_surfaces; + set m_reservoirs; + set m_flowDevices; + set m_walls; map m_counts; //!< Map used for default name generation unique_ptr m_integ; @@ -407,9 +412,6 @@ class ReactorNet : public FuncEval bool m_integrator_init = false; //!< True if integrator initialization is current size_t m_nv = 0; - //! m_start[n] is the starting point in the state vector for reactor n - vector m_start; - vector m_atol; double m_rtol = 1.0e-9; double m_rtolsens = 1.0e-4; diff --git a/include/cantera/zeroD/ReactorSurface.h b/include/cantera/zeroD/ReactorSurface.h index 9f11f1c7c1..85362fc111 100644 --- a/include/cantera/zeroD/ReactorSurface.h +++ b/include/cantera/zeroD/ReactorSurface.h @@ -7,6 +7,7 @@ #define CT_REACTOR_SURFACE_H #include "cantera/zeroD/ReactorBase.h" +#include "cantera/kinetics/InterfaceKinetics.h" namespace Cantera { @@ -54,7 +55,7 @@ class ReactorSurface : public ReactorBase } //! Accessor for the InterfaceKinetics object - Kinetics* kinetics() { + InterfaceKinetics* kinetics() { return m_kinetics.get(); } @@ -90,38 +91,162 @@ class ReactorSurface : public ReactorBase //! number of surface species. void getCoverages(double* cov) const; - //! Set the coverages and temperature in the surface phase object to the - //! values for this surface. The temperature is set to match the bulk phase - //! of the attached Reactor. - //! @since Prior to %Cantera 3.2, this operation was performed by syncState() - void restoreState() override; - - //! Set the coverages for this ReactorSurface based on the attached SurfPhase. - //! @since Behavior changed in %Cantera 3.2 for consistency with - //! ReactorBase::syncState(). Previously, this method performed the inverse - //! operation of setting the ReactorSurface state based on the SurfPhase and - //! attached Reactor. - void syncState() override; + void getState(double* y) override; + void initialize(double t0=0.0) override; + void updateState(double* y) override; + void eval(double t, double* LHS, double* RHS) override; void addSensitivityReaction(size_t rxn) override; + // void addSensitivitySpeciesEnthalpy(size_t k) override; + void applySensitivity(double* params) override; + void resetSensitivity(double* params) override; - //! Set reaction rate multipliers. `params` is the global vector of - //! sensitivity parameters. This function is called within - //! ReactorNet::eval() before the reaction rates are evaluated. - void setSensitivityParameters(const double* params); + size_t componentIndex(const string& nm) const override; + string componentName(size_t k) override; - //! Set reaction rate multipliers back to their initial values. This - //! function is called within ReactorNet::eval() after all rates have been - //! evaluated. - void resetSensitivityParameters(); + // vector steadyConstraints() const override; + // double upperBound(size_t k) const override; + // double lowerBound(size_t k) const override; + // void resetBadValues(double* y) override; + // void setDerivativeSettings(AnyMap& settings) override; protected: double m_area = 1.0; shared_ptr m_surf; - shared_ptr m_kinetics; + shared_ptr m_kinetics; vector m_reactors; - vector m_cov; +}; + +//! A surface where the state variables are the total number of moles of each species. +//! +//! This class provides the approximate Jacobian elements for interactions between +//! itself and the IdealGasMoleReactor and ConstPressureIdealGasMoleReactor classes +//! needed to work with the AdaptivePreconditioner class. +//! +//! @ingroup reactorGroup +//! @since New in %Cantera 4.0. +class MoleReactorSurface : public ReactorSurface +{ +public: + using ReactorSurface::ReactorSurface; + string type() const override { return "MoleReactorSurface"; } + void initialize(double t0=0.0) override; + void getState(double* y) override; + void updateState(double* y) override; + void eval(double t, double* LHS, double* RHS) override; + void getJacobianElements(vector>& trips) override; + +protected: + //! Temporary storage for coverages + vector m_cov_tmp; + + //! Temporary storage for d(moles)/d(moles) scaling factor + vector m_f_species; + + //! Temporary storage for d(energy)/d(moles) scaling factors + vector m_f_energy; + + //! Mapping from InterfaceKinetics species index to ReactorNet state vector index. + //! Set to -1 for species not included in the reactor network state vector. + vector::StorageIndex> m_kin2net; + + //! Mapping from InterfaceKinetics species index to the corresponding reactor. + //! Set to `nullptr` for surface species. + vector m_kin2reactor; + +}; + +//! A surface in contact with a FlowReactor. +//! +//! May represent the reactor wall or a catalytic surface within the reactor. +//! @ingroup reactorGroup +class FlowReactorSurface : public ReactorSurface +{ +public: + //! @copydoc ReactorSurface::ReactorSurface + FlowReactorSurface(shared_ptr soln, + const vector>& reactors, + bool clone, + const string& name="(none)"); + + string type() const override { + return "FlowReactorSurface"; + } + + void evalDae(double t, double* y, double* ydot, double* residual) override; + void getStateDae(double* y, double* ydot) override; + void getConstraints(double* constraints) override; + + //! Surface area per unit length [m] + //! @note If unspecified by the user, this will be calculated assuming the surface + //! is the wall of a cylindrical reactor. + double area() const override; + + //! Set the reactor surface area per unit length [m]. + //! + //! If the surface is the wall of the reactor, then this is the perimeter of the + //! reactor. If the surface represents something like a catalyst monolith, this is + //! the inverse of the surface area to volume ratio. + void setArea(double A) override { m_area = A; } + + //! Get the steady state tolerances used to determine the initial state for + //! surface coverages + double initialAtol() const { + return m_ss_atol; + } + + //! Set the steady state tolerances used to determine the initial state for + //! surface coverages + void setInitialAtol(double atol) { + m_ss_atol = atol; + } + + //! Get the steady state tolerances used to determine the initial state for + //! surface coverages + double initialRtol() const { + return m_ss_rtol; + } + + //! Set the steady state tolerances used to determine the initial state for + //! surface coverages + void setInitialRtol(double rtol) { + m_ss_rtol = rtol; + } + + //! Get the steady state tolerances used to determine the initial state for + //! surface coverages + double initialMaxSteps() const { + return m_max_ss_steps; + } + + //! Set the steady state tolerances used to determine the initial state for + //! surface coverages + void setInitialMaxSteps(int max_steps) { + m_max_ss_steps = max_steps; + } + + //! Get the steady state tolerances used to determine the initial state for + //! surface coverages + double initialMaxErrorFailures() const { + return m_max_ss_error_fails; + } + + //! Set the steady state tolerances used to determine the initial state for + //! surface coverages + void setInitialMaxErrorFailures(int max_fails) { + m_max_ss_error_fails = max_fails; + } + +protected: + //! steady-state relative tolerance, used to determine initial surface coverages + double m_ss_rtol = 1e-7; + //! steady-state absolute tolerance, used to determine initial surface coverages + double m_ss_atol = 1e-14; + //! maximum number of steady-state coverage integrator-steps + int m_max_ss_steps = 20000; + //! maximum number of steady-state integrator error test failures + int m_max_ss_error_fails = 10; }; } diff --git a/interfaces/cython/cantera/ctcxx.pxd b/interfaces/cython/cantera/ctcxx.pxd index 2083ab8cfa..7cb2c10626 100644 --- a/interfaces/cython/cantera/ctcxx.pxd +++ b/interfaces/cython/cantera/ctcxx.pxd @@ -21,6 +21,95 @@ ctypedef stdmap[string,double] Composition import numpy as np cimport numpy as np +from libcpp.vector cimport vector + +# See https://github.com/cython/cython/pull/6539 +# TODO: Replace once Cython >= 3.1.0 is required +cdef extern from "" namespace "std" nogil: + # Only Extent = std::dynamic_extent is supported until Cython can also + # support integer templates. See https://github.com/cython/cython/pull/426 + cdef cppclass span[T]: + ctypedef T value_type + ctypedef size_t size_type + ctypedef ptrdiff_t difference_type + + size_t extent + + cppclass iterator: + iterator() except + + iterator(iterator&) except + + T& operator*() + iterator operator++() + iterator operator--() + iterator operator++(int) + iterator operator--(int) + iterator operator+(size_type) + iterator operator-(size_type) + difference_type operator-(iterator) + difference_type operator-(const_iterator) + bint operator==(iterator) + bint operator==(const_iterator) + bint operator!=(iterator) + bint operator!=(const_iterator) + bint operator<(iterator) + bint operator<(const_iterator) + bint operator>(iterator) + bint operator>(const_iterator) + bint operator<=(iterator) + bint operator<=(const_iterator) + bint operator>=(iterator) + bint operator>=(const_iterator) + + cppclass reverse_iterator: + reverse_iterator() except + + reverse_iterator(reverse_iterator&) except + + T& operator*() + reverse_iterator operator++() + reverse_iterator operator--() + reverse_iterator operator++(int) + reverse_iterator operator--(int) + reverse_iterator operator+(size_type) + reverse_iterator operator-(size_type) + difference_type operator-(iterator) + difference_type operator-(const_iterator) + bint operator==(reverse_iterator) + bint operator==(const_reverse_iterator) + bint operator!=(reverse_iterator) + bint operator!=(const_reverse_iterator) + bint operator<(reverse_iterator) + bint operator<(const_reverse_iterator) + bint operator>(reverse_iterator) + bint operator>(const_reverse_iterator) + bint operator<=(reverse_iterator) + bint operator<=(const_reverse_iterator) + bint operator>=(reverse_iterator) + bint operator>=(const_reverse_iterator) + + span() + span(T*, size_type) except + # span[It](It, size_type) + span(T*, T*) except + # span[It, End](It, End) + span(vector&) # span[U, N](array[T, N]& arr) + span(span&) + + T& operator[](ssize_t) + + T& back() + iterator begin() + T* data() + bint empty() + iterator end() + span[T] first(size_type) + T& front() + span[T] last(size_type) + reverse_iterator rbegin() + reverse_iterator rend() + size_type size() + span[T] subspan(size_type) + span[T] subspan(size_type, size_type) + + cdef size_t dynamic_extent + + cdef extern from "cantera/cython/funcWrapper.h": cdef int translate_exception() diff --git a/interfaces/cython/cantera/delegator.pyx b/interfaces/cython/cantera/delegator.pyx index c5b01bfd53..19537b8762 100644 --- a/interfaces/cython/cantera/delegator.pyx +++ b/interfaces/cython/cantera/delegator.pyx @@ -274,7 +274,9 @@ cdef int assign_delegates(obj, CxxDelegator* delegator) except -1: cdef string cxx_name cdef string cxx_when obj._delegates = [] + valid_names = set() for name, options in obj.delegatable_methods.items(): + valid_names.add(name) if len(options) == 3: # Delegate with pre-selected mode, without using prefix on method name when = options[2] @@ -375,6 +377,14 @@ cdef int assign_delegates(obj, CxxDelegator* delegator) except -1: # collector obj._delegates.append(method) + # Check for methods on the base object that have no matching delegate + for name in dir(obj): + if (name.startswith("before_") or name.startswith("after_") or + name.startswith("replace_")): + base_name = name.split("_", 1)[1] + if base_name not in valid_names: + raise ValueError(f"'{base_name}' is not a delegatable method") + return 0 diff --git a/interfaces/cython/cantera/jacobians.pxd b/interfaces/cython/cantera/jacobians.pxd index a5c6c45b18..54786ff2c2 100644 --- a/interfaces/cython/cantera/jacobians.pxd +++ b/interfaces/cython/cantera/jacobians.pxd @@ -44,20 +44,17 @@ cdef extern from "cantera/numerics/SystemJacobianFactory.h" namespace "Cantera": cdef class SystemJacobian: @staticmethod cdef wrap(shared_ptr[CxxSystemJacobian]) - cdef set_cxx_object(self) cdef shared_ptr[CxxSystemJacobian] _base + cdef CxxSystemJacobian* jac cdef class EigenSparseJacobian(SystemJacobian): - cdef set_cxx_object(self) - cdef CxxEigenSparseJacobian* sparse_jac + pass cdef class EigenSparseDirectJacobian(EigenSparseJacobian): pass cdef class AdaptivePreconditioner(EigenSparseJacobian): - cdef set_cxx_object(self) - cdef CxxAdaptivePreconditioner* adaptive + pass cdef class BandedJacobian(SystemJacobian): - cdef set_cxx_object(self) - cdef CxxMultiJac* band_jac + pass diff --git a/interfaces/cython/cantera/jacobians.pyx b/interfaces/cython/cantera/jacobians.pyx index bc951af06f..d7f806ee6a 100644 --- a/interfaces/cython/cantera/jacobians.pyx +++ b/interfaces/cython/cantera/jacobians.pyx @@ -22,7 +22,7 @@ cdef class SystemJacobian: def _cinit(self, *args, **kwargs): self._base = newSystemJacobian(stringify(self._type)) - self.set_cxx_object() + self.jac = self._base.get() @staticmethod cdef wrap(shared_ptr[CxxSystemJacobian] base): @@ -45,12 +45,9 @@ cdef class SystemJacobian: cdef SystemJacobian jac jac = cls(init=False) jac._base = base - jac.set_cxx_object() + jac.jac = base.get() return jac - cdef set_cxx_object(self): - pass - property side: """ Get/Set the side of the system matrix where the preconditioner is applied. @@ -58,11 +55,10 @@ cdef class SystemJacobian: by all solver types. """ def __get__(self): - return pystr(self._base.get().preconditionerSide()) + return pystr(self.jac.preconditionerSide()) def __set__(self, side): - self._base.get().setPreconditionerSide(stringify(side)) - + self.jac.setPreconditionerSide(stringify(side)) cdef class EigenSparseJacobian(SystemJacobian): """ @@ -72,18 +68,15 @@ cdef class EigenSparseJacobian(SystemJacobian): _type = "eigen-sparse" - cdef set_cxx_object(self): - self.sparse_jac = self._base.get() - def print_contents(self): - self.sparse_jac.printPreconditioner() + (self.jac).printPreconditioner() property matrix: """ Property to retrieve the latest internal preconditioner matrix. """ def __get__(self): - cdef CxxSparseMatrix smat = self.sparse_jac.matrix() + cdef CxxSparseMatrix smat = (self.jac).matrix() return get_from_sparse(smat, smat.rows(), smat.cols()) property jacobian: @@ -91,7 +84,7 @@ cdef class EigenSparseJacobian(SystemJacobian): Property to retrieve the latest Jacobian. """ def __get__(self): - cdef CxxSparseMatrix smat = self.sparse_jac.jacobian() + cdef CxxSparseMatrix smat = (self.jac).jacobian() return get_from_sparse(smat, smat.rows(), smat.cols()) @@ -108,9 +101,6 @@ cdef class AdaptivePreconditioner(EigenSparseJacobian): _type = "Adaptive" linear_solver_type = "GMRES" - cdef set_cxx_object(self): - self.adaptive = self._base.get() - property threshold: """ The threshold of the preconditioner is used to remove or prune any off diagonal @@ -127,10 +117,10 @@ cdef class AdaptivePreconditioner(EigenSparseJacobian): Default is 0.0. """ def __get__(self): - return self.adaptive.threshold() + return (self.jac).threshold() def __set__(self, val): - self.adaptive.setThreshold(val) + (self.jac).setThreshold(val) property ilut_fill_factor: """ @@ -147,10 +137,10 @@ cdef class AdaptivePreconditioner(EigenSparseJacobian): Default is the state size divided by 4. """ def __set__(self, val): - self.adaptive.setIlutFillFactor(val) + (self.jac).setIlutFillFactor(val) def __get__(self): - return self.adaptive.ilutFillFactor() + return (self.jac).ilutFillFactor() property ilut_drop_tol: """ @@ -165,11 +155,10 @@ cdef class AdaptivePreconditioner(EigenSparseJacobian): Default is 1e-10. """ def __set__(self, val): - self.adaptive.setIlutDropTol(val) + (self.jac).setIlutDropTol(val) def __get__(self): - return self.adaptive.ilutDropTol() - + return (self.jac).ilutDropTol() cdef class BandedJacobian(SystemJacobian): """ @@ -178,6 +167,3 @@ cdef class BandedJacobian(SystemJacobian): """ _type = "banded-direct" linear_solver_type = "direct" - - cdef set_cxx_object(self): - self.band_jac = self._base.get() diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index 2db654e336..db8c9ea3fd 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -38,11 +38,14 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": CxxReactorBase() except +translate_exception string type() shared_ptr[CxxSolution] phase() - void restoreState() except +translate_exception void syncState() except +translate_exception double volume() string name() - void setName(string) + void setName(string) except +translate_exception + size_t neq() + size_t componentIndex(string&) except +translate_exception + string componentName(size_t) except +translate_exception + void getState(double*) except +translate_exception void setInitialVolume(double) void addSensitivityReaction(size_t) except +translate_exception @@ -52,10 +55,6 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": cbool chemistryEnabled() void setEnergyEnabled(cbool) cbool energyEnabled() - size_t componentIndex(string&) except +translate_exception - string componentName(size_t) except +translate_exception - size_t neq() - void getState(double*) except +translate_exception CxxSparseMatrix jacobian() except +translate_exception CxxSparseMatrix finiteDifferenceJacobian() except +translate_exception void setAdvanceLimit(string&, double) except +translate_exception @@ -67,29 +66,32 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": double speed() void setArea(double) except +translate_exception double area() except +translate_exception - void setSurfaceAreaToVolumeRatio(double) except +translate_exception - double surfaceAreaToVolumeRatio() except +translate_exception - double inletSurfaceAtol() - void setInletSurfaceAtol(double) except +translate_exception - double inletSurfaceRtol() - void setInletSurfaceRtol(double) except +translate_exception - double inletSurfaceMaxSteps() - void setInletSurfaceMaxSteps(int) except +translate_exception - int inletSurfaceMaxErrorFailures() - void setInletSurfaceMaxErrorFailures(int) except +translate_exception + # reactor surface - cdef cppclass CxxReactorSurface "Cantera::ReactorSurface": - string type() - string name() - void setName(string) except +translate_exception + cdef cppclass CxxReactorSurface "Cantera::ReactorSurface" (CxxReactorBase): double area() void setArea(double) void setKinetics(CxxKinetics*) except +translate_exception void setCoverages(double*) void setCoverages(Composition&) except +translate_exception - void syncState() + + cdef cppclass CxxFlowReactorSurface "Cantera::FlowReactorSurface" (CxxReactorSurface): + CxxFlowReactorSurface() except +translate_exception + + double initialAtol() + void setInitialAtol(double) except +translate_exception + double initialRtol() + void setInitialRtol(double) except +translate_exception + double initialMaxSteps() + void setInitialMaxSteps(int) except +translate_exception + int initialMaxErrorFailures() + void setInitialMaxErrorFailures(int) except +translate_exception + + cdef shared_ptr[CxxReactorBase] CxxNewReactorSurface "newReactorSurface" ( + string, shared_ptr[CxxSolution], vector[shared_ptr[CxxReactorBase]]&, + cbool, string) except +translate_exception cdef shared_ptr[CxxReactorBase] CxxNewReactorSurface "newReactorSurface" ( shared_ptr[CxxSolution], vector[shared_ptr[CxxReactorBase]]&, cbool, string) except +translate_exception @@ -214,8 +216,7 @@ cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera": void setExpansionRate(double) double heatRate() void setHeatRate(double) - void restoreThermoState() except +translate_exception - void restoreSurfaceState(size_t) except +translate_exception + span[double] surfaceProductionRates() ctypedef CxxReactorAccessor* CxxReactorAccessorPtr @@ -272,11 +273,28 @@ cdef class FlowReactor(Reactor): cdef class ExtensibleReactor(Reactor): cdef public _delegates cdef CxxReactorAccessor* accessor + cdef public double[:] surface_production_rates + """ + An array containing the total production rates [kmol/s] of bulk species on + surfaces. Updated from adjacent surfaces as part of ``Reactor.eval``. + """ cdef class ReactorSurface(ReactorBase): cdef CxxReactorSurface* surface cdef list _reactors +cdef class FlowReactorSurface(ReactorSurface): + pass + +cdef class ExtensibleReactorSurface(ReactorSurface): + cdef public _delegates + cdef public double[:] surface_production_rates + """ + An array containing the net production rates [kmol/m²/s] of surface and bulk species + on the surface, with the order corresponding to the `InterfaceKinetics` object. + Updated as part of ``ReactorSurface.update_state``. + """ + cdef class ConnectorNode: cdef shared_ptr[CxxConnectorNode] _node cdef CxxConnectorNode* node diff --git a/interfaces/cython/cantera/reactor.pyi b/interfaces/cython/cantera/reactor.pyi index 96d43a0222..e175da66f3 100644 --- a/interfaces/cython/cantera/reactor.pyi +++ b/interfaces/cython/cantera/reactor.pyi @@ -38,6 +38,11 @@ class ReactorBase: def name(self) -> str: ... @name.setter def name(self, name: str) -> None: ... + def component_index(self, name: str) -> int: ... + def component_name(self, i: int) -> str: ... + @property + def n_vars(self) -> int: ... + def get_state(self) -> Array: ... def syncState(self) -> None: ... @property def phase(self) -> Solution: ... @@ -99,11 +104,6 @@ class Reactor(ReactorBase): @energy_enabled.setter def energy_enabled(self, value: bool) -> None: ... def add_sensitivity_species_enthalpy(self, k: int) -> None: ... - def component_index(self, name: str) -> int: ... - def component_name(self, i: int) -> str: ... - @property - def n_vars(self) -> int: ... - def get_state(self) -> Array: ... @property def jacobian(self) -> Array: ... @property @@ -129,26 +129,6 @@ class FlowReactor(Reactor): @area.setter def area(self, area: float) -> None: ... @property - def inlet_surface_atol(self) -> float: ... - @inlet_surface_atol.setter - def inlet_surface_atol(self, atol: float) -> None: ... - @property - def inlet_surface_rtol(self) -> float: ... - @inlet_surface_rtol.setter - def inlet_surface_rtol(self, rtol: float) -> None: ... - @property - def inlet_surface_max_steps(self) -> int: ... - @inlet_surface_max_steps.setter - def inlet_surface_max_steps(self, nsteps: int) -> None: ... - @property - def inlet_surface_max_error_failures(self) -> int: ... - @inlet_surface_max_error_failures.setter - def inlet_surface_max_error_failures(self, nsteps: int) -> None: ... - @property - def surface_area_to_volume_ratio(self) -> float: ... - @surface_area_to_volume_ratio.setter - def surface_area_to_volume_ratio(self, sa_to_vol: float) -> None: ... - @property def speed(self) -> float: ... class ExtensibleReactor(Reactor): @@ -167,8 +147,7 @@ class ExtensibleReactor(Reactor): def heat_rate(self) -> float: ... @heat_rate.setter def heat_rate(self, qdot: float) -> None: ... - def restore_thermo_state(self) -> None: ... - def restore_surface_state(self, n: int) -> None: ... + surface_production_rates: Array class ExtensibleIdealGasReactor(ExtensibleReactor): ... class ExtensibleConstPressureReactor(ExtensibleReactor): ... @@ -204,6 +183,7 @@ class ReactorSurface: def reactor(self) -> Reactor: ... @property def reactors(self) -> list[Reactor]: ... + surface_production_rates: Array def draw( self, graph: Digraph | None = None, @@ -216,6 +196,31 @@ class ReactorSurface: species_units: Literal["percent", "ppm"] = "percent", ) -> Digraph: ... +class FlowReactorSurface(ReactorSurface): + @property + def initial_atol(self) -> float: ... + @initial_atol.setter + def initial_atol(self, atol: float) -> None: ... + @property + def initial_rtol(self) -> float: ... + @initial_rtol.setter + def initial_rtol(self, rtol: float) -> None: ... + @property + def initial_max_steps(self) -> int: ... + @initial_max_steps.setter + def initial_max_steps(self, max_steps: int) -> None: ... + @property + def initial_max_error_failures(self) -> int: ... + @initial_max_error_failures.setter + def initial_max_error_failures(self, n: int) -> None: ... + +class ExtensibleReactorSurface(ReactorSurface): + delegatable_methods: dict[str, tuple[str, str]] + @property + def n_vars(self) -> int: ... + +class ExtensibleMoleReactorSurface(ExtensibleReactorSurface): ... + class ConnectorNode: node_type: ClassVar[str] edge_attr: dict[str, str] diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 921cbeebb0..45e7b5e635 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -17,7 +17,7 @@ cdef class ReactorBase: reactor_type = "none" def __cinit__(self, _SolutionBase phase, *args, name="(none)", clone=True, **kwargs): - if self.reactor_type != "ReactorSurface": + if not isinstance(self, ReactorSurface): self._rbase = newReactorBase(stringify(self.reactor_type), phase._base, clone, stringify(name)) self.rbase = self._rbase.get() @@ -47,11 +47,79 @@ cdef class ReactorBase: def __set__(self, name): self.rbase.setName(stringify(name)) + def component_index(self, name): + """ + Returns the index of the component named ``name`` in the system. This determines + the index of the component in the vector of sensitivity coefficients. ``name`` + is either a species name or the name of a reactor state variable, for example + ``'int_energy'`` or ``'temperature'``, depending on the reactor's equations. + """ + k = self.rbase.componentIndex(stringify(name)) + if k == -1: + raise IndexError('No such component: {!r}'.format(name)) + return k + + def component_name(self, int i): + """ + Returns the name of the component with index ``i`` within the array of + variables returned by `get_state`. This is the inverse of + `component_index`. + """ + return pystr(self.rbase.componentName(i)) + + property n_vars: + """ + The number of state variables in the reactor. + Equal to: + + `Reactor` and `IdealGasReactor`: `n_species` + 3 (mass, volume, + internal energy or temperature). + + `ConstPressureReactor` and `IdealGasConstPressureReactor`: + `n_species` + 2 (mass, enthalpy or temperature). + """ + def __get__(self): + return self.rbase.neq() + + def get_state(self): + """ + Get the state vector of the reactor. + + The order of the variables (that is, rows) is: + + `Reactor` or `IdealGasReactor`: + + - 0 - mass + - 1 - volume + - 2 - internal energy or temperature + - 3+ - mass fractions of the species + + `ConstPressureReactor` or `IdealGasConstPressureReactor`: + + - 0 - mass + - 1 - enthalpy or temperature + - 2+ - mass fractions of the species + + You can use the function `component_index` to determine the location of + a specific component from its name, or `component_name` to determine the + name from the index. + """ + if not self.n_vars: + raise CanteraError('Reactor empty or network not initialized.') + cdef np.ndarray[np.double_t, ndim=1] y = np.zeros(self.n_vars) + self.rbase.getState(&y[0]) + return y + def syncState(self): """ Set the state of the Reactor to match that of the associated `ThermoPhase` object. After calling syncState(), call ReactorNet.reinitialize() before further integration. + + .. deprecated:: 4.0 + Manual synchronization of reactor state is no longer required. Call + `ReactorNet.reinitialize` directly to indicate a change in state that + requires integrator reinitialization. """ self.rbase.syncState() @@ -63,7 +131,6 @@ cdef class ReactorBase: .. versionchanged:: 3.2 Renamed from ``thermo``. """ - self.rbase.restoreState() return self._phase property volume: @@ -287,69 +354,6 @@ cdef class Reactor(ReactorBase): """ self.reactor.addSensitivitySpeciesEnthalpy(self.phase.species_index(k)) - def component_index(self, name): - """ - Returns the index of the component named ``name`` in the system. This determines - the index of the component in the vector of sensitivity coefficients. ``name`` - is either a species name or the name of a reactor state variable, for example - ``'int_energy'`` or ``'temperature'``, depending on the reactor's equations. - """ - k = self.reactor.componentIndex(stringify(name)) - if k == -1: - raise IndexError('No such component: {!r}'.format(name)) - return k - - def component_name(self, int i): - """ - Returns the name of the component with index ``i`` within the array of - variables returned by `get_state`. This is the inverse of - `component_index`. - """ - return pystr(self.reactor.componentName(i)) - - property n_vars: - """ - The number of state variables in the reactor. - Equal to: - - `Reactor` and `IdealGasReactor`: `n_species` + 3 (mass, volume, - internal energy or temperature). - - `ConstPressureReactor` and `IdealGasConstPressureReactor`: - `n_species` + 2 (mass, enthalpy or temperature). - """ - def __get__(self): - return self.reactor.neq() - - def get_state(self): - """ - Get the state vector of the reactor. - - The order of the variables (that is, rows) is: - - `Reactor` or `IdealGasReactor`: - - - 0 - mass - - 1 - volume - - 2 - internal energy or temperature - - 3+ - mass fractions of the species - - `ConstPressureReactor` or `IdealGasConstPressureReactor`: - - - 0 - mass - - 1 - enthalpy or temperature - - 2+ - mass fractions of the species - - You can use the function `component_index` to determine the location of - a specific component from its name, or `component_name` to determine the - name from the index. - """ - if not self.n_vars: - raise CanteraError('Reactor empty or network not initialized.') - cdef np.ndarray[np.double_t, ndim=1] y = np.zeros(self.n_vars) - self.reactor.getState(&y[0]) - return y - property jacobian: """ Get the local, reactor-specific Jacobian or an approximation thereof @@ -486,63 +490,6 @@ cdef class FlowReactor(Reactor): def area(self, area): (self.reactor).setArea(area) - @property - def inlet_surface_atol(self): - """ - Get/Set the steady-state tolerances used to determine the initial surface - species coverages. - """ - return (self.reactor).inletSurfaceAtol() - - @inlet_surface_atol.setter - def inlet_surface_atol(self, atol): - (self.reactor).setInletSurfaceAtol(atol) - - @property - def inlet_surface_rtol(self): - """ - Get/Set the steady-state tolerances used to determine the initial surface - species coverages. - """ - return (self.reactor).inletSurfaceRtol() - - @inlet_surface_rtol.setter - def inlet_surface_rtol(self, rtol): - (self.reactor).setInletSurfaceRtol(rtol) - - @property - def inlet_surface_max_steps(self): - """ - Get/Set the maximum number of integrator steps used to determine the initial - surface species coverages. - """ - return (self.reactor).inletSurfaceMaxSteps() - - @inlet_surface_max_steps.setter - def inlet_surface_max_steps(self, nsteps): - (self.reactor).setInletSurfaceMaxSteps(nsteps) - - @property - def inlet_surface_max_error_failures(self): - """ - Get/Set the maximum number of integrator error failures allowed when determining - the initial surface species coverages. - """ - return (self.reactor).inletSurfaceMaxErrorFailures() - - @inlet_surface_max_error_failures.setter - def inlet_surface_max_error_failures(self, nsteps): - (self.reactor).setInletSurfaceMaxErrorFailures(nsteps) - - @property - def surface_area_to_volume_ratio(self): - """Get/Set the surface area to volume ratio of the reactor [1/m]""" - return (self.reactor).surfaceAreaToVolumeRatio() - - @surface_area_to_volume_ratio.setter - def surface_area_to_volume_ratio(self, sa_to_vol): - (self.reactor).setSurfaceAreaToVolumeRatio(sa_to_vol) - @property def speed(self): """ Speed [m/s] of the flow in the reactor at the current position """ @@ -554,6 +501,10 @@ cdef class ExtensibleReactor(Reactor): A base class for a reactor with delegated methods where the base functionality corresponds to the `Reactor` class. + The ``__init__`` method of the derived class should allocate and size any + internal variables and set the total number of state variables associated with this + reactor, `n_vars` (if it is different from the base class). + The following methods of the C++ :ct:`Reactor` class can be modified by a Python class which inherits from this class. For each method, the name below should be prefixed with ``before_``, ``after_``, or ``replace_``, indicating @@ -568,16 +519,11 @@ cdef class ExtensibleReactor(Reactor): from the supplied method and the base class method. ``initialize(self, t0: double) -> None`` - Responsible for allocating and setting the sizes of any internal - variables, initializing attached walls, and setting the total number of - state variables associated with this reactor, `n_vars`. + Responsible for initialization that can only be performed after connecting the + elements of the reactor network, such as initializing attached walls. Called once before the start of time integration. - ``sync_state(self) -> None`` - Responsible for setting the state of the reactor to correspond to the - state of the associated ThermoPhase object. - ``get_state(self, y : double[:]) -> None`` Responsible for populating the state vector ``y`` (length `n_vars`) with the initial state of the reactor. @@ -586,16 +532,6 @@ cdef class ExtensibleReactor(Reactor): Responsible for setting the state of the reactor object from the values in the state vector ``y`` (length `n_vars`) - ``update_surface_state(self, y : double[:]) -> None`` - Responsible for setting the state of surface phases in this reactor - from the values in the state vector ``y``. The length of ``y`` is the - total number of surface species in all surfaces. - - ``get_surface_initial_conditions(self, y : double[:]) -> None`` - Responsible for populating the state vector ``y`` with the initial - state of each surface phase in this reactor. The length of ``y`` is the - total number of surface species in all surfaces. - ``update_connected(self, update_pressure : bool) -> None`` Responsible for storing properties which may be accessed by connected reactors, and for updating the mass flow rates of connected flow devices. @@ -612,44 +548,31 @@ cdef class ExtensibleReactor(Reactor): 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`` - Responsible for calculating the ``LHS`` and ``RHS`` (length: total number of - surface species in all surfaces) of the ODEs for surface species coverages, - and the molar production rate of bulk phase species ``sdot`` (length: number - of bulk phase species). - ``component_name(i : int) -> string`` Returns the name of the state vector component with index ``i`` ``component_index(name: string) -> int`` Returns the index of the state vector component named ``name`` - - ``species_index(name : string) -> int`` - Returns the index of the species named ``name``, in either the bulk - phase or a surface phase, relative to the start of the species terms in - the state vector. """ reactor_type = "ExtensibleReactor" delegatable_methods = { 'initialize': ('initialize', 'void(double)'), - 'sync_state': ('syncState', 'void()'), 'get_state': ('getState', 'void(double*)'), 'update_state': ('updateState', 'void(double*)'), - 'update_surface_state': ('updateSurfaceState', 'void(double*)'), - 'get_surface_initial_conditions': ('getSurfaceInitialConditions', 'void(double*)'), 'update_connected': ('updateConnected', 'void(bool)'), 'eval': ('eval', 'void(double, double*, double*)'), 'eval_walls': ('evalWalls', 'void(double)'), - 'eval_surfaces': ('evalSurfaces', 'void(double*,double*,double*)'), 'component_name': ('componentName', 'string(size_t)'), 'component_index': ('componentIndex', 'size_t(string)'), - 'species_index': ('speciesIndex', 'size_t(string)') } def __cinit__(self, *args, **kwargs): self.accessor = dynamic_cast[CxxReactorAccessorPtr](self.rbase) + cdef span[double] sdot = \ + dynamic_cast[CxxReactorAccessorPtr](self.rbase).surfaceProductionRates() + self.surface_production_rates = sdot.data() def __init__(self, *args, **kwargs): assign_delegates(self, dynamic_cast[CxxDelegatorPtr](self.rbase)) @@ -690,20 +613,6 @@ cdef class ExtensibleReactor(Reactor): def heat_rate(self, qdot): self.accessor.setHeatRate(qdot) - def restore_thermo_state(self): - """ - Set the state of the thermo object to correspond to the state of the - reactor. - """ - self.accessor.restoreThermoState() - - def restore_surface_state(self, n): - """ - Set the state of the thermo object for surface ``n`` to correspond to the - state of that surface - """ - self.accessor.restoreSurfaceState(n) - cdef class ExtensibleIdealGasReactor(ExtensibleReactor): """ @@ -806,7 +715,8 @@ cdef class ReactorSurface(ReactorBase): """ reactor_type = "ReactorSurface" - def __cinit__(self, _SolutionBase phase, r, clone=True, name="(none)", **kwargs): + def __cinit__(self, _SolutionBase phase, r, *, clone=True, name="(none)", + kind=None, **kwargs): cdef ReactorBase adj cdef vector[shared_ptr[CxxReactorBase]] cxx_adj if isinstance(r, ReactorBase): @@ -824,11 +734,20 @@ cdef class ReactorSurface(ReactorBase): raise TypeError("Parameter 'r' should be a ReactorBase object or a list " "of ReactorBase objects.") - self._rbase = CxxNewReactorSurface(phase._base, cxx_adj, clone, stringify(name)) + if kind is None and self.reactor_type != "ReactorSurface": + kind = self.reactor_type + + if kind is not None: + self._rbase = CxxNewReactorSurface(stringify(kind), phase._base, cxx_adj, + clone, stringify(name)) + else: + self._rbase = CxxNewReactorSurface(phase._base, cxx_adj, clone, + stringify(name)) + self.rbase = self._rbase.get() self.surface = (self.rbase) - def __init__(self, phase=None, r=None, *, clone=True, + def __init__(self, phase=None, r=None, *, kind=None, clone=True, name="(none)", A=None, node_attr=None): super().__init__(phase, name=name) @@ -850,7 +769,6 @@ cdef class ReactorSurface(ReactorBase): def __get__(self): if self._phase is None: raise CanteraError('No kinetics manager present') - self.rbase.restoreState() return self._phase.coverages def __set__(self, coverages): if self._phase is None: @@ -933,6 +851,109 @@ cdef class ReactorSurface(ReactorBase): print_state, species, species_units) +cdef class FlowReactorSurface(ReactorSurface): + reactor_type = "FlowReactorSurface" + + @property + def initial_atol(self): + """ + Get/Set the steady-state tolerances used to determine the initial surface + species coverages. + """ + return (self.surface).initialAtol() + + @initial_atol.setter + def initial_atol(self, atol): + (self.surface).setInitialAtol(atol) + @property + def initial_rtol(self): + """ + Get/Set the steady-state tolerances used to determine the initial surface + species coverages. + """ + return (self.surface).initialRtol() + + @initial_rtol.setter + def initial_rtol(self, rtol): + (self.surface).setInitialRtol(rtol) + + @property + def initial_max_steps(self): + """ + Get/Set the maximum number of integrator steps used to determine the initial + surface species coverages. + """ + return (self.surface).initialMaxSteps() + + @initial_max_steps.setter + def initial_max_steps(self, nsteps): + (self.surface).setInitialMaxSteps(nsteps) + + @property + def initial_max_error_failures(self): + """ + Get/Set the maximum number of integrator error failures allowed when determining + the initial surface species coverages. + """ + return (self.surface).initialMaxErrorFailures() + + @initial_max_error_failures.setter + def initial_max_error_failures(self, nsteps): + (self.surface).setInitialMaxErrorFailures(nsteps) + + +cdef class ExtensibleReactorSurface(ReactorSurface): + """ + A base class for a reactor surface with delegated methods where the base + functionality corresponds to the `ReactorSurface` class. + + Methods of the base class can be modified in the same way as in class + `ExtensibleReactor`. The methods that can be modified are:: + + - ``initialize(self, t0: double) -> None`` + - ``get_state(self, y : double[:]) -> None`` + - ``update_state(self, y : double[:]) -> None`` + - ``eval(self, t : double, LHS : double[:], RHS : double[:]) -> None`` + - ``component_name(i : int) -> string`` + - ``component_index(name: string) -> int`` + """ + + reactor_type = "ExtensibleReactorSurface" + + delegatable_methods = { + 'initialize': ('initialize', 'void(double)'), + 'get_state': ('getState', 'void(double*)'), + 'update_state': ('updateState', 'void(double*)'), + 'eval': ('eval', 'void(double, double*, double*)'), + 'component_name': ('componentName', 'string(size_t)'), + 'component_index': ('componentIndex', 'size_t(string)'), + } + + def __init__(self, *args, **kwargs): + assign_delegates(self, dynamic_cast[CxxDelegatorPtr](self.rbase)) + cdef span[double] sdot = \ + dynamic_cast[CxxReactorAccessorPtr](self.rbase).surfaceProductionRates() + self.surface_production_rates = sdot.data() + super().__init__(*args, **kwargs) + + property n_vars: + """ + Get/Set the number of state variables in the reactor. + """ + def __get__(self): + return self.reactor.neq() + def __set__(self, n): + dynamic_cast[CxxReactorAccessorPtr](self.rbase).setNEq(n) + + +cdef class ExtensibleMoleReactorSurface(ExtensibleReactorSurface): + """ + A variant of `ExtensibleReactorSurface` where the base behavior corresponds to the + :ct:`MoleReactorSurface` class. + """ + reactor_type = "ExtensibleMoleReactorSurface" + + cdef class ConnectorNode: """ Common base class for walls and flow devices. diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctreactor.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctreactor.yaml index 666c39be8f..f1d7e4f25c 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctreactor.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctreactor.yaml @@ -14,7 +14,9 @@ recipes: - name: new wraps: newReactorBase - name: newSurface - wraps: newReactorSurface + wraps: newReactorSurface(shared_ptr, const vector>&, bool, const string&) +- name: newSurfaceOfType + wraps: newReactorSurface(const string&, shared_ptr, const vector>&, bool, const string&) - name: type - name: name - name: setName diff --git a/samples/cxx/openmp_ignition/openmp_ignition.cpp b/samples/cxx/openmp_ignition/openmp_ignition.cpp index 82b748e876..aaeb1f4650 100644 --- a/samples/cxx/openmp_ignition/openmp_ignition.cpp +++ b/samples/cxx/openmp_ignition/openmp_ignition.cpp @@ -66,7 +66,6 @@ void run() // Set up the problem gas->setState_TPX(T0[i], OneAtm, "CH4:0.5, O2:1.0, N2:3.76"); - reactor.syncState(); net.setInitialTime(0.0); // Integrate until we satisfy a crude estimate of the ignition delay diff --git a/samples/python/reactors/custom2.py b/samples/python/reactors/custom2.py index 1215139b8a..7094ffe360 100644 --- a/samples/python/reactors/custom2.py +++ b/samples/python/reactors/custom2.py @@ -31,8 +31,7 @@ def __init__(self, *args, neighbor, **kwargs): self.k_wall = 1e-2 # proportionality constant, a_wall = k_wall * delta P self.neighbor = neighbor - def after_initialize(self, t0): - # The initialize function for the base Reactor class will have set + # The constructor for the base Reactor class will have set # n_vars to already include the volume, internal energy, mass, and mass # fractions of all the species. Increase this by one to account for # the added variable of the wall velocity. diff --git a/samples/python/reactors/pfr.py b/samples/python/reactors/pfr.py index f96173ca76..92a7009dbf 100644 --- a/samples/python/reactors/pfr.py +++ b/samples/python/reactors/pfr.py @@ -127,7 +127,6 @@ for n in range(n_steps): # Set the state of the reservoir to match that of the previous reactor upstream.phase.TDY = r2.phase.TDY - upstream.syncState() # integrate the reactor forward in time until steady state is reached sim2.reinitialize() sim2.advance_to_steady_state() diff --git a/samples/python/reactors/porous_media_burner.py b/samples/python/reactors/porous_media_burner.py index 35a9bdb054..13a7aec723 100644 --- a/samples/python/reactors/porous_media_burner.py +++ b/samples/python/reactors/porous_media_burner.py @@ -195,8 +195,6 @@ def __init__(self, *args, props, **kwargs): self.V = self.A * props.length # reactor volume (m^3) self.molecular_weights = None self.species_offset = None - - def after_initialize(self, t0): self.n_vars += 1 # additional equation for the solid temperature self.index_Ts = self.n_vars - 1 self.species_offset = self.component_index(self.phase.species_name(0)) diff --git a/samples/python/reactors/surf_pfr.py b/samples/python/reactors/surf_pfr.py index 9835f42c3d..1769f4c541 100644 --- a/samples/python/reactors/surf_pfr.py +++ b/samples/python/reactors/surf_pfr.py @@ -52,14 +52,16 @@ class and the SUNDIALS IDA solver, in contrast to the approximation as a chain o # create a new reactor r = ct.FlowReactor(gas) r.area = area -r.surface_area_to_volume_ratio = cat_area_per_vol * porosity r.mass_flow_rate = mass_flow_rate r.energy_enabled = False # Add the reacting surface to the reactor -rsurf = ct.ReactorSurface(surf, r) +rsurf = ct.FlowReactorSurface(surf, r) +#r.surface_area_to_volume_ratio = cat_area_per_vol * porosity +rsurf.area = cat_area_per_vol * porosity * area sim = ct.ReactorNet([r]) +sim.verbose = True output_data = [] n = 0 diff --git a/samples/python/reactors/surf_pfr_chain.py b/samples/python/reactors/surf_pfr_chain.py index cebf2bee2f..33be1765a9 100644 --- a/samples/python/reactors/surf_pfr_chain.py +++ b/samples/python/reactors/surf_pfr_chain.py @@ -17,6 +17,7 @@ import csv import cantera as ct +ct.CanteraError.set_stack_trace_depth(10) # unit conversion factors to SI cm = 0.01 @@ -105,8 +106,7 @@ for n in range(NReactors): # Set the state of the reservoir to match that of the previous reactor - gas.TDY = r.phase.TDY - upstream.syncState() + upstream.phase.TDY = r.phase.TDY sim.reinitialize() sim.advance_to_steady_state() dist = n * rlen * 1.0e3 # distance in mm diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 4b8ca10108..276bc288f5 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -652,8 +652,8 @@ void InterfaceKinetics::setDerivativeSettings(const AnyMap& settings) void InterfaceKinetics::getDerivativeSettings(AnyMap& settings) const { - settings["skip-coverage-dependence"] = m_jac_skip_electrochemistry; - settings["skip-electrochemistry"] = m_jac_skip_coverage_dependence; + settings["skip-coverage-dependence"] = m_jac_skip_coverage_dependence; + settings["skip-electrochemistry"] = m_jac_skip_electrochemistry; settings["rtol-delta"] = m_jac_rtol_delta; } diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index e7f7046e80..f266d26484 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -319,6 +319,17 @@ size_t Kinetics::kineticsSpeciesIndex(const string& nm, bool raise) const return npos; } +size_t Kinetics::speciesOffset(const ThermoPhase& phase) const +{ + for (size_t n = 0; n < m_thermo.size(); n++) { + if (&phase == m_thermo[n].get()) { + return m_start[n]; + } + } + throw CanteraError("Kinetics::speciesOffset", + "Phase '{}' not found in Kinetics object", phase.name()); +} + ThermoPhase& Kinetics::speciesPhase(const string& nm) { for (size_t n = 0; n < m_thermo.size(); n++) { diff --git a/src/zeroD/ConstPressureMoleReactor.cpp b/src/zeroD/ConstPressureMoleReactor.cpp index ea3aa745d9..134eb9263f 100644 --- a/src/zeroD/ConstPressureMoleReactor.cpp +++ b/src/zeroD/ConstPressureMoleReactor.cpp @@ -15,27 +15,27 @@ namespace Cantera { +ConstPressureMoleReactor::ConstPressureMoleReactor(shared_ptr sol, + const string& name) + : ConstPressureMoleReactor(sol, true, name) +{ +} + +ConstPressureMoleReactor::ConstPressureMoleReactor(shared_ptr sol, bool clone, + const string& name) + : MoleReactor(sol, clone, name) +{ + m_nv = 1 + m_nsp; // enthalpy and moles of each species +} + void ConstPressureMoleReactor::getState(double* y) { - if (m_thermo == 0) { - throw CanteraError("ConstPressureMoleReactor::getState", - "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); // set mass to be used in getMoles function m_mass = m_thermo->density() * m_vol; // set the first array element to enthalpy y[0] = m_thermo->enthalpy_mass() * m_thermo->density() * m_vol; // get moles of species in remaining state getMoles(y + m_sidx); - // set the remaining components to the surface species moles on the walls - getSurfaceInitialConditions(y+m_nsp+m_sidx); -} - -void ConstPressureMoleReactor::initialize(double t0) -{ - MoleReactor::initialize(t0); - m_nv -= 1; // const pressure system loses 1 more variable from MoleReactor } void ConstPressureMoleReactor::updateState(double* y) @@ -52,7 +52,6 @@ void ConstPressureMoleReactor::updateState(double* y) } m_vol = m_mass / m_thermo->density(); updateConnected(false); - updateSurfaceState(y + m_nsp + m_sidx); } void ConstPressureMoleReactor::eval(double time, double* LHS, double* RHS) @@ -60,8 +59,7 @@ void ConstPressureMoleReactor::eval(double time, double* LHS, double* RHS) double* dndt = RHS + m_sidx; // kmol per s evalWalls(time); - - m_thermo->restoreState(m_state); + updateSurfaceProductionRates(); const vector& imw = m_thermo->inverseMolecularWeights(); @@ -69,9 +67,6 @@ void ConstPressureMoleReactor::eval(double time, double* LHS, double* RHS) m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" } - // evaluate reactor surfaces - evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data()); - // external heat transfer double dHdt = m_Qdot; @@ -113,7 +108,7 @@ size_t ConstPressureMoleReactor::componentIndex(const string& nm) const return 0; } try { - return speciesIndex(nm) + m_sidx; + return m_thermo->speciesIndex(nm) + m_sidx; } catch (const CanteraError&) { throw CanteraError("ConstPressureMoleReactor::componentIndex", "Component '{}' not found", nm); @@ -124,22 +119,11 @@ string ConstPressureMoleReactor::componentName(size_t k) { if (k == 0) { return "enthalpy"; } else if (k >= m_sidx && k < neq()) { - k -= m_sidx; - if (k < m_thermo->nSpecies()) { - return m_thermo->speciesName(k); - } else { - k -= m_thermo->nSpecies(); - } - for (auto& S : m_surfaces) { - ThermoPhase* th = S->thermo(); - if (k < th->nSpecies()) { - return th->speciesName(k); - } else { - k -= th->nSpecies(); - } - } + return m_thermo->speciesName(k - m_sidx); + } else { + throw IndexError("ConstPressureMoleReactor::componentName", + "component", k, m_nv); } - throw IndexError("ConstPressureMoleReactor::componentName", "component", k, m_nv); } double ConstPressureMoleReactor::upperBound(size_t k) const { diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index b479218c3e..0d2eedcf18 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -14,14 +14,21 @@ namespace Cantera { -void ConstPressureReactor::getState(double* y) +ConstPressureReactor::ConstPressureReactor(shared_ptr sol, + const string& name) + : ConstPressureReactor(sol, true, name) { - if (m_thermo == 0) { - throw CanteraError("ConstPressureReactor::getState", - "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); +} +ConstPressureReactor::ConstPressureReactor(shared_ptr sol, bool clone, + const string& name) + : Reactor(sol, clone, name) +{ + m_nv = 2 + m_nsp; // mass, enthalpy, and mass fractions of each species +} + +void ConstPressureReactor::getState(double* y) +{ // set the first component to the total mass y[0] = m_thermo->density() * m_vol; @@ -31,15 +38,6 @@ void ConstPressureReactor::getState(double* y) // set components y+2 ... y+K+1 to the mass fractions Y_k of each species m_thermo->getMassFractions(y+2); - // set the remaining components to the surface species - // coverages on the walls - getSurfaceInitialConditions(y + m_nsp + 2); -} - -void ConstPressureReactor::initialize(double t0) -{ - Reactor::initialize(t0); - m_nv -= 1; // Constant pressure reactor has one fewer state variable } void ConstPressureReactor::updateState(double* y) @@ -56,7 +54,6 @@ void ConstPressureReactor::updateState(double* y) } m_vol = m_mass / m_thermo->density(); updateConnected(false); - updateSurfaceState(y + m_nsp + 2); } void ConstPressureReactor::eval(double time, double* LHS, double* RHS) @@ -67,12 +64,10 @@ void ConstPressureReactor::eval(double time, double* LHS, double* RHS) dmdt = 0.0; evalWalls(time); + updateSurfaceProductionRates(); - m_thermo->restoreState(m_state); const vector& mw = m_thermo->molecularWeights(); const double* Y = m_thermo->massFractions(); - - evalSurfaces(LHS + m_nsp + 2, RHS + m_nsp + 2, m_sdot.data()); double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin()); dmdt += mdot_surf; @@ -137,7 +132,7 @@ size_t ConstPressureReactor::componentIndex(const string& nm) const return 1; } try { - return speciesIndex(nm) + 2; + return m_thermo->speciesIndex(nm) + 2; } catch (const CanteraError&) { throw CanteraError("ConstPressureReactor::componentIndex", "Component '{}' not found", nm); @@ -150,20 +145,7 @@ string ConstPressureReactor::componentName(size_t k) { } else if (k == 1) { return "enthalpy"; } else if (k >= 2 && k < neq()) { - k -= 2; - if (k < m_thermo->nSpecies()) { - return m_thermo->speciesName(k); - } else { - k -= m_thermo->nSpecies(); - } - for (auto& S : m_surfaces) { - ThermoPhase* th = S->thermo(); - if (k < th->nSpecies()) { - return th->speciesName(k); - } else { - k -= th->nSpecies(); - } - } + return m_thermo->speciesName(k - 2); } throw IndexError("ConstPressureReactor::componentName", "component", k, m_nv); } diff --git a/src/zeroD/FlowReactor.cpp b/src/zeroD/FlowReactor.cpp index 045713aff8..7667f9c8a0 100644 --- a/src/zeroD/FlowReactor.cpp +++ b/src/zeroD/FlowReactor.cpp @@ -17,17 +17,28 @@ namespace Cantera { +FlowReactor::FlowReactor(shared_ptr sol, const string& name) + : FlowReactor(sol, true, name) +{ +} + +FlowReactor::FlowReactor(shared_ptr sol, bool clone, const string& name) + : IdealGasReactor(sol, clone, name) +{ + m_nv = 4 + m_nsp; // rho, u, P, T, and species mass fractions + m_rho = m_thermo->density(); + // resize temporary arrays + m_wdot.resize(m_nsp); + m_hk.resize(m_nsp); +} + void FlowReactor::getStateDae(double* y, double* ydot) { - if (m_thermo == nullptr) { - throw CanteraError("FlowReactor::getStateDae", "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); m_thermo->getMassFractions(y+m_offset_Y); const vector& mw = m_thermo->molecularWeights(); // set the first component to the initial density - y[0] = m_rho; + y[0] = m_thermo->density(); if (m_u < 0) { throw CanteraError("FlowReactor::getStateDae", @@ -38,35 +49,20 @@ void FlowReactor::getStateDae(double* y, double* ydot) y[1] = m_u; // set the third component to the initial pressure - y[2] = m_P; + y[2] = m_thermo->pressure(); // set the fourth component to the initial temperature - y[3] = m_T; + y[3] = m_thermo->temperature(); if (m_chem) { m_kin->getNetProductionRates(m_wdot.data()); // "omega dot" } - // need to advance the reactor surface to steady state to get the initial - // coverages - for (auto m_surf : m_surfaces) { - m_surf->restoreState(); - auto kin = static_cast(m_surf->kinetics()); - kin->advanceCoverages(100.0, m_ss_rtol, m_ss_atol, 0, m_max_ss_steps, - m_max_ss_error_fails); - auto& surf = dynamic_cast(kin->thermo(0)); - vector cov(surf.nSpecies()); - surf.getCoverages(cov.data()); - m_surf->setCoverages(cov.data()); - } - // set the initial coverages - getSurfaceInitialConditions(y + m_offset_Y + m_nsp); + updateSurfaceProductionRates(); // reset ydot vector std::fill(ydot, ydot + m_nv, 0.0); - // calculate initial d(coverage)/dt values - evalSurfaces(ydot + m_nsp + 4, m_sdot.data()); // next, we must solve for the initial derivative values // a . ydot = b @@ -80,11 +76,11 @@ void FlowReactor::getStateDae(double* y, double* ydot) DenseMatrix a(m_offset_Y + m_nsp, m_offset_Y + m_nsp); // first row is the ideal gas equation, differentiated with respect to z - a(0, 0) = - GasConstant * m_T / m_thermo->meanMolecularWeight(); + a(0, 0) = - GasConstant * m_thermo->temperature() / m_thermo->meanMolecularWeight(); a(0, 2) = 1; a(0, 3) = - m_rho * GasConstant / m_thermo->meanMolecularWeight(); for (size_t i = 0; i < m_nsp; ++i) { - a(0, m_offset_Y + i) = - m_rho * m_T / mw[i]; + a(0, m_offset_Y + i) = - m_rho * m_thermo->temperature() / mw[i]; } // initialize the second row from mass conservation equation (Kee 16.48) @@ -109,9 +105,8 @@ void FlowReactor::getStateDae(double* y, double* ydot) // get (perim / Ac) * sum(sk' * wk), used multiple places double h_sk_wk = 0; - double hydraulic = surfaceAreaToVolumeRatio(); for (size_t i = 0; i < m_nsp; ++i) { - h_sk_wk += hydraulic * m_sdot[i] * mw[i]; + h_sk_wk += m_sdot[i] * mw[i] / m_area; } // RHS of ideal gas eq. is zero @@ -135,7 +130,7 @@ void FlowReactor::getStateDae(double* y, double* ydot) // h_mass * Wk = h_mol m_thermo->getPartialMolarEnthalpies(m_hk.data()); for (size_t i = 0; i < m_nsp; ++i) { - ydot[3] -= m_hk[i] * (m_wdot[i] + hydraulic * m_sdot[i]); + ydot[3] -= m_hk[i] * (m_wdot[i] + m_sdot[i] / m_area); } } else { ydot[3] = 0; @@ -145,66 +140,20 @@ void FlowReactor::getStateDae(double* y, double* ydot) // - Yk * (perim / Ac) * sum(sk' * Wk) + wk' * Wk + (perim / Ac) * sk' * Wk for (size_t i = 0; i < m_nsp; ++i) { ydot[m_offset_Y + i] = -y[m_offset_Y + i] * h_sk_wk + - mw[i] * (m_wdot[i] + hydraulic * m_sdot[i]); + mw[i] * (m_wdot[i] + m_sdot[i] / m_area); } // and solve solve(a, ydot, 1, 0); } -void FlowReactor::initialize(double t0) -{ - Reactor::initialize(t0); - m_thermo->restoreState(m_state); - // initialize state - m_T = m_thermo->temperature(); - m_rho = m_thermo->density(); - m_P = m_thermo->pressure(); - m_T = m_thermo->temperature(); - // resize temporary arrays - m_wdot.resize(m_nsp); - m_hk.resize(m_nsp); - // set number of variables to the number of non-species equations - // i.e., density, velocity, pressure and temperature - // plus the number of species in the gas phase - m_nv = m_offset_Y + m_nsp; - if (m_surfaces.size()) { - size_t n_surf_species = 0; - size_t n_total_species = 0; - for (auto const &m_surf : m_surfaces) { - // finally, add the number of species the surface for the site-fraction - // equations - n_surf_species += m_surf->thermo()->nSpecies(); - n_total_species += m_surf->kinetics()->nTotalSpecies(); - } - m_nv += n_surf_species; - m_sdot_temp.resize(n_total_species); - } -} - -void FlowReactor::syncState() -{ - Reactor::syncState(); - m_rho = m_thermo->density(); - m_P = m_thermo->pressure(); - m_T = m_thermo->temperature(); -} - void FlowReactor::updateState(double* y) { // Set the mass fractions and density of the mixture. m_rho = y[0]; m_u = y[1]; - m_P = y[2]; - m_T = y[3]; - double* mss = y + m_offset_Y; - m_thermo->setMassFractions_NoNorm(mss); - m_thermo->setState_TP(m_T, m_P); - - // update surface - updateSurfaceState(y + m_nsp + m_offset_Y); - - m_thermo->saveState(m_state); + m_thermo->setMassFractions_NoNorm(y + m_offset_Y); + m_thermo->setState_TP(y[3], y[2]); } void FlowReactor::setMassFlowRate(double mdot) @@ -219,39 +168,13 @@ void FlowReactor::setArea(double area) { setMassFlowRate(mdot); } -double FlowReactor::surfaceAreaToVolumeRatio() const { - if (m_sa_to_vol > 0) { - return m_sa_to_vol; - } - - // assuming a cylinder, volume = Pi * r^2 * L, and perimeter = 2 * Pi * r * L - // where L is the length, and r is the radius of the reactor - // hence, perimeter / area = 2 * Pi * r * L / (Pi * L * r^2) = 2 / r - return 2.0 / sqrt(m_area / Pi); -} - -void FlowReactor::updateSurfaceState(double* y) -{ - size_t loc = 0; - for (auto& S : m_surfaces) { - // set the coverages without normalization to avoid over-constraining - // the system. - // note: the ReactorSurface class doesn't normalize when calling setCoverages - S->setCoverages(y+loc); - S->restoreState(); - loc += S->thermo()->nSpecies(); - } -} - void FlowReactor::evalDae(double time, double* y, double* ydot, double* residual) { - m_thermo->restoreState(m_state); - - evalSurfaces(ydot + m_nsp + 4, m_sdot.data()); + updateSurfaceProductionRates(); const vector& mw = m_thermo->molecularWeights(); double sk_wk = 0; for (size_t i = 0; i < m_nsp; ++i) { - sk_wk = m_sdot[i] * mw[i]; + sk_wk += m_sdot[i] * mw[i] / m_area; } m_thermo->getPartialMolarEnthalpies(m_hk.data()); // get net production @@ -270,13 +193,12 @@ void FlowReactor::evalDae(double time, double* y, double* ydot, double* residual //! use mass continuity for velocity residual //! Kee.'s Chemically Reacting Flow, Eq. 16.48 - double hydraulic = surfaceAreaToVolumeRatio(); - residual[1] = m_u * drhodz + m_rho * dudz - sk_wk * hydraulic; + residual[1] = m_u * drhodz + m_rho * dudz - sk_wk; //! Use conservation of momentum for pressure residual //! Kee.'s Chemically Reacting Flow, Eq. 16.62 //! [NOTE: friction is neglected] - residual[2] = m_rho * m_u * dudz + m_u * hydraulic * sk_wk + dPdz; + residual[2] = m_rho * m_u * dudz + m_u * sk_wk + dPdz; //! use conservation of energy for temperature residual //! Kee.'s Chemically Reacting Flow, Eq. 16.58 @@ -289,7 +211,7 @@ void FlowReactor::evalDae(double time, double* y, double* ydot, double* residual double cp_mass = m_thermo->cp_mass(); residual[3] = m_rho * m_u * cp_mass * dTdz; for (size_t i = 0; i < m_nsp; ++i) { - residual[3] += m_hk[i] * (m_wdot[i] + hydraulic * m_sdot[i]); + residual[3] += m_hk[i] * (m_wdot[i] + m_sdot[i] / m_area); } } else { residual[3] = dTdz; @@ -300,49 +222,23 @@ void FlowReactor::evalDae(double time, double* y, double* ydot, double* residual double dSumYdz = 0; for (size_t i = 0; i < m_nsp; ++i) { residual[i + m_offset_Y] = m_rho * m_u * ydot[i + m_offset_Y] + - y[i + m_offset_Y] * hydraulic * sk_wk - - mw[i] * (m_wdot[i] + hydraulic * m_sdot[i]); + y[i + m_offset_Y] * sk_wk - + mw[i] * (m_wdot[i] + m_sdot[i] / m_area); dSumYdz += ydot[i + m_offset_Y]; } // Spread d/dz(sum(Y)) = 0 constraint across all species equations. `scale` is // defined to make the size of the error in sum(Y) comparable to the overall rtol. - double scale = 0.1 * m_rho * m_u / m_ss_rtol; + // TODO: find a better way to access a value of rtol. + double rtol = 1e-7; + double scale = 0.1 * m_rho * m_u / rtol; for (size_t i = 0; i < m_nsp; ++i) { residual[i + m_offset_Y] += scale * std::max(0.0, y[i + m_offset_Y]) * dSumYdz; } - - // surface algebraic constraints - if (m_surfaces.size()) { - size_t loc = m_offset_Y + m_nsp; // offset into residual vector - for (auto &m_surf : m_surfaces) { - Kinetics* kin = m_surf->kinetics(); - size_t nk = m_surf->thermo()->nSpecies(); - kin->getNetProductionRates(m_sdot_temp.data()); - double sum = y[loc]; - for (size_t i = 1; i < nk; ++i) { - //! net surface production rate residuals - //! Kee.'s Chemically Reacting Flow, Eq. 16.63 - residual[loc + i] = m_sdot_temp[i]; - //! next, evaluate the algebraic constraint to explicitly conserve - //! surface site fractions. - //! Kee.'s Chemically Reacting Flow, Eq. 16.64 - sum += y[loc + i]; - } - // set residual of the first species sum - 1 to ensure - // surface site conservation. - // Note: the first species is selected to be consistent with - // Reactor::evalSurfaces - residual[loc] = sum - 1.0; - loc += nk; - } - } } void FlowReactor::getConstraints(double* constraints) { // mark all variables differential equations unless otherwise specified std::fill(constraints, constraints + m_nv, 1.0); - // the species coverages are algebraic constraints - std::fill(constraints + m_offset_Y + m_nsp, constraints + m_nv, 0.0); } size_t FlowReactor::componentIndex(const string& nm) const @@ -357,7 +253,7 @@ size_t FlowReactor::componentIndex(const string& nm) const return 3; } try { - return speciesIndex(nm) + m_offset_Y; + return m_thermo->speciesIndex(nm) + m_offset_Y; } catch (const CanteraError&) { throw CanteraError("FlowReactor::componentIndex", "Component '{}' not found", nm); @@ -375,20 +271,7 @@ string FlowReactor::componentName(size_t k) } else if (k == 3) { return "temperature"; } else if (k >= 4 && k < neq()) { - k -= 4; - if (k < m_thermo->nSpecies()) { - return m_thermo->speciesName(k); - } else { - k -= m_thermo->nSpecies(); - } - for (auto& S : m_surfaces) { - ThermoPhase* th = S->thermo(); - if (k < th->nSpecies()) { - return th->speciesName(k); - } else { - k -= th->nSpecies(); - } - } + return m_thermo->speciesName(k - 4); } throw IndexError("FlowReactor::componentName", "component", k, m_nv); } diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index 648dea2748..564ad8ce9a 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -19,19 +19,12 @@ namespace Cantera void IdealGasConstPressureMoleReactor::getState(double* y) { - if (m_thermo == 0) { - throw CanteraError("IdealGasConstPressureMoleReactor::getState", - "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); // get mass for calculations m_mass = m_thermo->density() * m_vol; // set the first component to the temperature y[0] = m_thermo->temperature(); // get moles of species in remaining state getMoles(y + m_sidx); - // set the remaining components to the surface species moles on the walls - getSurfaceInitialConditions(y + m_nsp + m_sidx); } void IdealGasConstPressureMoleReactor::initialize(double t0) @@ -53,8 +46,9 @@ void IdealGasConstPressureMoleReactor::updateState(double* y) m_thermo->setMolesNoTruncate(y + m_sidx); m_thermo->setState_TP(y[0], m_pressure); m_vol = m_mass / m_thermo->density(); + m_thermo->getPartialMolarEnthalpies(m_hk.data()); + m_TotalCp = m_mass * m_thermo->cp_mass(); updateConnected(false); - updateSurfaceState(y + m_nsp + m_sidx); } void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RHS) @@ -63,19 +57,13 @@ void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RH double* dndt = RHS + m_sidx; // kmol per s evalWalls(time); - - m_thermo->restoreState(m_state); - - m_thermo->getPartialMolarEnthalpies(&m_hk[0]); + updateSurfaceProductionRates(); const vector& imw = m_thermo->inverseMolecularWeights(); if (m_chem) { m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" } - // evaluate reactor surfaces - evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data()); - // external heat transfer mcpdTdt += m_Qdot; @@ -108,71 +96,30 @@ void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RH } if (m_energy) { - LHS[0] = m_mass * m_thermo->cp_mass(); + LHS[0] = m_TotalCp; } else { RHS[0] = 0.0; } } -Eigen::SparseMatrix IdealGasConstPressureMoleReactor::jacobian() +void IdealGasConstPressureMoleReactor::getJacobianElements( + vector>& trips) { - if (m_nv == 0) { - throw CanteraError("IdealGasConstPressureMoleReactor::jacobian", - "Reactor must be initialized first."); - } - // clear former jacobian elements - m_jac_trips.clear(); // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as // d (dot(omega)) / d c_j, it is later transformed appropriately. Eigen::SparseMatrix dnk_dnj = m_kin->netProductionRates_ddCi(); - // species size that accounts for surface species - size_t ssize = m_nv - m_sidx; - // map derivatives from the surface chemistry jacobian - // to the reactor jacobian - if (!m_surfaces.empty()) { - vector> species_trips(dnk_dnj.nonZeros()); - for (int k = 0; k < dnk_dnj.outerSize(); k++) { - for (Eigen::SparseMatrix::InnerIterator it(dnk_dnj, k); it; ++it) { - species_trips.emplace_back(static_cast(it.row()), - static_cast(it.col()), it.value()); - } - } - addSurfaceJacobian(species_trips); - dnk_dnj.resize(ssize, ssize); - dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end()); - } - // get net production rates - Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize); - // gas phase net production rates - m_kin->getNetProductionRates(netProductionRates.data()); - // surface phase net production rates mapped to reactor gas phase - for (auto &S: m_surfaces) { - auto curr_kin = S->kinetics(); - vector prod_rates(curr_kin->nTotalSpecies()); - curr_kin->getNetProductionRates(prod_rates.data()); - for (size_t i = 0; i < curr_kin->nTotalSpecies(); i++) { - try { - size_t row = speciesIndex(curr_kin->kineticsSpeciesName(i)); - netProductionRates[row] += prod_rates[i]; - } catch (...) { - // species do not map - } - } - } - double molarVol = m_thermo->molarVolume(); + // add species to species derivatives elements to the jacobian // calculate ROP derivatives, excluding the terms -n_i / (V * N) dc_i/dn_j // as it substantially reduces matrix sparsity + size_t offset = static_cast(m_offset + m_sidx); + // double molarVol = m_thermo->molarVolume(); for (int k = 0; k < dnk_dnj.outerSize(); k++) { for (Eigen::SparseMatrix::InnerIterator it(dnk_dnj, k); it; ++it) { - // gas phase species need the addition of V / N * omega_dot - if (static_cast(it.row()) < m_nsp) { - it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol; - } - m_jac_trips.emplace_back(static_cast(it.row() + m_sidx), - static_cast(it.col() + m_sidx), it.value()); + trips.emplace_back(it.row() + offset, it.col() + offset, it.value()); } } + // Temperature Derivatives if (m_energy) { // getting perturbed state for finite difference @@ -198,39 +145,40 @@ Eigen::SparseMatrix IdealGasConstPressureMoleReactor::jacobian() for (size_t j = 0; j < m_nv; j++) { double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j]; double ydotCurrent = rhsCurrent[j] / lhsCurrent[j]; - m_jac_trips.emplace_back(static_cast(j), 0, + trips.emplace_back(static_cast(j + m_offset), static_cast(m_offset), (ydotPerturbed - ydotCurrent) / deltaTemp); } // d T_dot/dnj // allocate vectors for whole system - Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize); - Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize); + Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(m_nsp); + Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(m_nsp); + Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(m_nsp); // gas phase m_thermo->getPartialMolarCp(specificHeat.data()); m_thermo->getPartialMolarEnthalpies(enthalpy.data()); + m_kin->getNetProductionRates(netProductionRates.data()); // scale production rates by the volume for gas species for (size_t i = 0; i < m_nsp; i++) { netProductionRates[i] *= m_vol; } double qdot = enthalpy.dot(netProductionRates); - // find denominator ahead of time - double NCp = 0.0; - double* moles = yCurrent.data() + m_sidx; - for (size_t i = 0; i < ssize; i++) { - NCp += moles[i] * specificHeat[i]; - } - double denom = 1 / (NCp * NCp); + double denom = 1 / (m_TotalCp * m_TotalCp); Eigen::VectorXd hk_dnkdnj_sums = dnk_dnj.transpose() * enthalpy; // Add derivatives to jac by spanning columns - for (size_t j = 0; j < ssize; j++) { - m_jac_trips.emplace_back(0, static_cast(j + m_sidx), - (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom); + for (size_t j = 0; j < m_nsp; j++) { + trips.emplace_back(m_offset, static_cast(j + m_offset + m_sidx), + (specificHeat[j] * qdot - m_TotalCp * hk_dnkdnj_sums[j]) * denom); } } - // convert triplets to sparse matrix - Eigen::SparseMatrix jac(m_nv, m_nv); - jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); - return jac; +} + +void IdealGasConstPressureMoleReactor::getJacobianScalingFactors( + double& f_species, double* f_energy) +{ + f_species = 1.0 / m_vol; + for (size_t k = 0; k < m_nsp; k++) { + f_energy[k] = - m_hk[k] / m_TotalCp; + } } size_t IdealGasConstPressureMoleReactor::componentIndex(const string& nm) const @@ -239,7 +187,7 @@ size_t IdealGasConstPressureMoleReactor::componentIndex(const string& nm) const return 0; } try { - return speciesIndex(nm) + m_sidx; + return m_thermo->speciesIndex(nm) + m_sidx; } catch (const CanteraError&) { throw CanteraError("IdealGasConstPressureReactor::componentIndex", "Component '{}' not found", nm); @@ -250,20 +198,7 @@ string IdealGasConstPressureMoleReactor::componentName(size_t k) { if (k == 0) { return "temperature"; } else if (k >= m_sidx && k < neq()) { - k -= m_sidx; - if (k < m_thermo->nSpecies()) { - return m_thermo->speciesName(k); - } else { - k -= m_thermo->nSpecies(); - } - for (auto& S : m_surfaces) { - ThermoPhase* th = S->thermo(); - if (k < th->nSpecies()) { - return th->speciesName(k); - } else { - k -= th->nSpecies(); - } - } + return m_thermo->speciesName(k - m_sidx); } throw IndexError("IdealGasConstPressureMoleReactor::componentName", "components", k, m_nv); diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 583479171a..d04a2a7163 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -15,12 +15,6 @@ namespace Cantera void IdealGasConstPressureReactor::getState(double* y) { - if (m_thermo == 0) { - throw CanteraError("IdealGasConstPressureReactor::getState", - "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); - // set the first component to the total mass y[0] = m_thermo->density() * m_vol; @@ -29,10 +23,6 @@ void IdealGasConstPressureReactor::getState(double* y) // set components y+2 ... y+K+1 to the mass fractions Y_k of each species m_thermo->getMassFractions(y+2); - - // set the remaining components to the surface species - // coverages on the walls - getSurfaceInitialConditions(y + m_nsp + 2); } void IdealGasConstPressureReactor::initialize(double t0) @@ -56,7 +46,6 @@ void IdealGasConstPressureReactor::updateState(double* y) m_thermo->setState_TP(y[1], m_pressure); m_vol = m_mass / m_thermo->density(); updateConnected(false); - updateSurfaceState(y + m_nsp + 2); } void IdealGasConstPressureReactor::eval(double time, double* LHS, double* RHS) @@ -69,15 +58,11 @@ void IdealGasConstPressureReactor::eval(double time, double* LHS, double* RHS) mcpdTdt = 0.0; evalWalls(time); - - m_thermo->restoreState(m_state); + updateSurfaceProductionRates(); const vector& mw = m_thermo->molecularWeights(); const double* Y = m_thermo->massFractions(); - - evalSurfaces(LHS + m_nsp + 2, RHS + m_nsp + 2, m_sdot.data()); double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin()); dmdt += mdot_surf; - m_thermo->getPartialMolarEnthalpies(&m_hk[0]); if (m_chem) { @@ -148,7 +133,7 @@ size_t IdealGasConstPressureReactor::componentIndex(const string& nm) const return 1; } try { - return speciesIndex(nm) + 2; + return m_thermo->speciesIndex(nm) + 2; } catch (const CanteraError&) { throw CanteraError("IdealGasConstPressureReactor::componentIndex", "Component '{}' not found", nm); diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index d7d39a7140..9c832210cf 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -19,12 +19,6 @@ namespace Cantera void IdealGasMoleReactor::getState(double* y) { - if (m_thermo == 0) { - throw CanteraError("IdealGasMoleReactor::getState", - "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); - // get mass for calculations m_mass = m_thermo->density() * m_vol; @@ -36,9 +30,6 @@ void IdealGasMoleReactor::getState(double* y) // get moles of species in remaining state getMoles(y + m_sidx); - // set the remaining components to the surface species moles on - // the walls - getSurfaceInitialConditions(y + m_nsp + m_sidx); } size_t IdealGasMoleReactor::componentIndex(const string& nm) const @@ -50,7 +41,7 @@ size_t IdealGasMoleReactor::componentIndex(const string& nm) const return 1; } try { - return speciesIndex(nm) + m_sidx; + return m_thermo->speciesIndex(nm) + m_sidx; } catch (const CanteraError&) { throw CanteraError("IdealGasMoleReactor::componentIndex", "Component '{}' not found", nm); @@ -119,8 +110,9 @@ void IdealGasMoleReactor::updateState(double* y) // set state m_thermo->setMolesNoTruncate(y + m_sidx); m_thermo->setState_TD(y[0], m_mass / m_vol); + m_thermo->getPartialMolarIntEnergies(m_uk.data()); + m_TotalCv = m_mass * m_thermo->cv_mass(); updateConnected(true); - updateSurfaceState(y + m_nsp + m_sidx); } void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS) @@ -129,19 +121,13 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS) double* dndt = RHS + m_sidx; // kmol per s evalWalls(time); - - m_thermo->restoreState(m_state); - - m_thermo->getPartialMolarIntEnergies(&m_uk[0]); + updateSurfaceProductionRates(); const vector& imw = m_thermo->inverseMolecularWeights(); if (m_chem) { m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" } - // evaluate surfaces - evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data()); - // external heat transfer and compression work mcvdTdt += - m_pressure * m_vdot + m_Qdot; @@ -177,48 +163,28 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS) RHS[1] = m_vdot; if (m_energy) { - LHS[0] = m_mass * m_thermo->cv_mass(); + LHS[0] = m_TotalCv; } else { RHS[0] = 0; } } -Eigen::SparseMatrix IdealGasMoleReactor::jacobian() +void IdealGasMoleReactor::getJacobianElements(vector>& trips) { - if (m_nv == 0) { - throw CanteraError("IdealGasMoleReactor::jacobian", - "Reactor must be initialized first."); - } - // clear former jacobian elements - m_jac_trips.clear(); // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as // d (dot(omega)) / d c_j, it is later transformed appropriately. Eigen::SparseMatrix dnk_dnj = m_kin->netProductionRates_ddCi(); - // species size that accounts for surface species - size_t ssize = m_nv - m_sidx; - // map derivatives from the surface chemistry jacobian - // to the reactor jacobian - if (!m_surfaces.empty()) { - vector> species_trips; - for (int k = 0; k < dnk_dnj.outerSize(); k++) { - for (Eigen::SparseMatrix::InnerIterator it(dnk_dnj, k); it; ++it) { - species_trips.emplace_back(static_cast(it.row()), - static_cast(it.col()), it.value()); - } - } - addSurfaceJacobian(species_trips); - dnk_dnj.resize(ssize, ssize); - dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end()); - } + // add species to species derivatives elements to the jacobian // calculate ROP derivatives, excluding the terms -n_i / (V * N) dc_i/dn_j // as it substantially reduces matrix sparsity for (int k = 0; k < dnk_dnj.outerSize(); k++) { for (Eigen::SparseMatrix::InnerIterator it(dnk_dnj, k); it; ++it) { - m_jac_trips.emplace_back(static_cast(it.row() + m_sidx), - static_cast(it.col() + m_sidx), it.value()); + trips.emplace_back(static_cast(it.row() + m_offset + m_sidx), + static_cast(it.col() + m_offset + m_sidx), it.value()); } } + // Temperature Derivatives if (m_energy) { // getting perturbed state for finite difference @@ -243,13 +209,13 @@ Eigen::SparseMatrix IdealGasMoleReactor::jacobian() for (size_t j = 0; j < m_nv; j++) { double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j]; double ydotCurrent = rhsCurrent[j] / lhsCurrent[j]; - m_jac_trips.emplace_back(static_cast(j), 0, - (ydotPerturbed - ydotCurrent) / deltaTemp); + trips.emplace_back(static_cast(j + m_offset), m_offset, + (ydotPerturbed - ydotCurrent) / deltaTemp); } // d T_dot/dnj - Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize); - Eigen::VectorXd internal_energy = Eigen::VectorXd::Zero(ssize); - Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize); + Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(m_nsp); + Eigen::VectorXd internal_energy = Eigen::VectorXd::Zero(m_nsp); + Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(m_nsp); // getting species data m_thermo->getPartialMolarIntEnergies(internal_energy.data()); m_kin->getNetProductionRates(netProductionRates.data()); @@ -261,25 +227,23 @@ Eigen::SparseMatrix IdealGasMoleReactor::jacobian() } // scale net production rates by volume to get molar rate double qdot = internal_energy.dot(netProductionRates); - // find the sum of n_i and cp_i - double NCv = 0.0; - double* moles = yCurrent.data() + m_sidx; - for (size_t i = 0; i < ssize; i++) { - NCv += moles[i] * specificHeat[i]; - } - // make denominator beforehand - double denom = 1 / (NCv * NCv); + double denom = 1 / (m_TotalCv * m_TotalCv); Eigen::VectorXd uk_dnkdnj_sums = dnk_dnj.transpose() * internal_energy; // add derivatives to jacobian - for (size_t j = 0; j < ssize; j++) { - m_jac_trips.emplace_back(0, static_cast(j + m_sidx), - (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom); + for (size_t j = 0; j < m_nsp; j++) { + trips.emplace_back(m_offset, static_cast(j + m_offset + m_sidx), + (specificHeat[j] * qdot - m_TotalCv * uk_dnkdnj_sums[j]) * denom); } } - // convert triplets to sparse matrix - Eigen::SparseMatrix jac(m_nv, m_nv); - jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); - return jac; +} + +void IdealGasMoleReactor::getJacobianScalingFactors( + double& f_species, double* f_energy) +{ + f_species = 1.0 / m_vol; + for (size_t k = 0; k < m_nsp; k++) { + f_energy[k] = - m_uk[k] / m_TotalCv; + } } } diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 890f5bb0f2..ddb8dc9883 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -15,12 +15,6 @@ namespace Cantera void IdealGasReactor::getState(double* y) { - if (m_thermo == 0) { - throw CanteraError("IdealGasReactor::getState", - "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); - // set the first component to the total mass m_mass = m_thermo->density() * m_vol; y[0] = m_mass; @@ -33,10 +27,6 @@ void IdealGasReactor::getState(double* y) // set components y+3 ... y+K+2 to the mass fractions of each species m_thermo->getMassFractions(y+3); - - // set the remaining components to the surface species - // coverages on the walls - getSurfaceInitialConditions(y + m_nsp + 3); } void IdealGasReactor::initialize(double t0) @@ -61,7 +51,6 @@ void IdealGasReactor::updateState(double* y) m_thermo->setMassFractions_NoNorm(y+3); m_thermo->setState_TD(y[2], m_mass / m_vol); updateConnected(true); - updateSurfaceState(y + m_nsp + 3); } void IdealGasReactor::eval(double time, double* LHS, double* RHS) @@ -71,7 +60,7 @@ void IdealGasReactor::eval(double time, double* LHS, double* RHS) double* mdYdt = RHS + 3; // mass * dY/dt evalWalls(time); - m_thermo->restoreState(m_state); + updateSurfaceProductionRates(); m_thermo->getPartialMolarIntEnergies(&m_uk[0]); const vector& mw = m_thermo->molecularWeights(); const double* Y = m_thermo->massFractions(); @@ -80,7 +69,6 @@ void IdealGasReactor::eval(double time, double* LHS, double* RHS) m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" } - evalSurfaces(LHS + m_nsp + 3, RHS + m_nsp + 3, m_sdot.data()); double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin()); dmdt += mdot_surf; @@ -157,7 +145,7 @@ size_t IdealGasReactor::componentIndex(const string& nm) const return 2; } try { - return speciesIndex(nm) + 3; + return m_thermo->speciesIndex(nm) + 3; } catch (const CanteraError&) { throw CanteraError("IdealGasReactor::componentIndex", "Component '{}' not found", nm); diff --git a/src/zeroD/MoleReactor.cpp b/src/zeroD/MoleReactor.cpp index ba10068b70..32c717832e 100644 --- a/src/zeroD/MoleReactor.cpp +++ b/src/zeroD/MoleReactor.cpp @@ -18,124 +18,15 @@ namespace bmt = boost::math::tools; namespace Cantera { -void MoleReactor::getSurfaceInitialConditions(double* y) +MoleReactor::MoleReactor(shared_ptr sol, const string& name) + : MoleReactor(sol, true, name) { - size_t loc = 0; - for (auto& S : m_surfaces) { - double area = S->area(); - auto currPhase = S->thermo(); - size_t tempLoc = currPhase->nSpecies(); - double surfDensity = currPhase->siteDensity(); - S->getCoverages(y + loc); - // convert coverages to moles - for (size_t i = 0; i < tempLoc; i++) { - y[i + loc] = y[i + loc] * area * surfDensity / currPhase->size(i); - } - loc += tempLoc; - } -} - -void MoleReactor::initialize(double t0) -{ - Reactor::initialize(t0); - m_nv -= 1; // moles gives the state one fewer variables -} - -void MoleReactor::updateSurfaceState(double* y) -{ - size_t loc = 0; - vector coverages(m_nv_surf, 0.0); - for (auto& S : m_surfaces) { - auto surf = S->thermo(); - double invArea = 1/S->area(); - double invSurfDensity = 1/surf->siteDensity(); - size_t tempLoc = surf->nSpecies(); - for (size_t i = 0; i < tempLoc; i++) { - coverages[i + loc] = y[i + loc] * invArea * surf->size(i) * invSurfDensity; - } - S->setCoverages(coverages.data()+loc); - loc += tempLoc; - } -} - -void MoleReactor::evalSurfaces(double* LHS, double* RHS, double* sdot) -{ - fill(sdot, sdot + m_nsp, 0.0); - size_t loc = 0; // offset into ydot - for (auto S : m_surfaces) { - Kinetics* kin = S->kinetics(); - SurfPhase* surf = S->thermo(); - double wallarea = S->area(); - size_t nk = surf->nSpecies(); - S->restoreState(); - kin->getNetProductionRates(&m_work[0]); - for (size_t k = 0; k < nk; k++) { - RHS[loc + k] = m_work[k] * wallarea / surf->size(k); - } - loc += nk; - - size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0)); - - for (size_t k = 0; k < m_nsp; k++) { - sdot[k] += m_work[bulkloc + k] * wallarea; - } - } } -void MoleReactor::addSurfaceJacobian(vector> &triplets) +MoleReactor::MoleReactor(shared_ptr sol, bool clone, const string& name) + : Reactor(sol, clone, name) { - size_t offset = m_nsp; - for (auto& S : m_surfaces) { - S->restoreState(); - double A = S->area(); - auto kin = S->kinetics(); - size_t nk = S->thermo()->nSpecies(); - // index of gas and surface phases to check if the species is in gas or surface - size_t spi = 0; - size_t gpi = kin->speciesPhaseIndex(kin->kineticsSpeciesIndex( - m_thermo->speciesName(0))); - // get surface jacobian in concentration units - Eigen::SparseMatrix surfJac = kin->netProductionRates_ddCi(); - // loop through surface specific jacobian and add elements to triplets vector - // accordingly - for (int k=0; k::InnerIterator it(surfJac, k); it; ++it) { - size_t row = it.row(); - size_t col = it.col(); - auto& rowPhase = kin->speciesPhase(row); - auto& colPhase = kin->speciesPhase(col); - size_t rpi = kin->phaseIndex(rowPhase.name(), true); - size_t cpi = kin->phaseIndex(colPhase.name(), true); - // check if the reactor kinetics object contains both phases to avoid - // any solid phases which may be included then use phases to map surf - // kinetics indicies to reactor kinetic indices - if ((rpi == spi || rpi == gpi) && (cpi == spi || cpi == gpi) ) { - // subtract start of phase - row -= kin->kineticsSpeciesIndex(0, rpi); - col -= kin->kineticsSpeciesIndex(0, cpi); - // since the gas phase is the first phase in the reactor state - // vector add the offset only if it is a surf species index to - // both row and col - row = (rpi != spi) ? row : row + offset; - // determine appropriate scalar to account for dimensionality - // gas phase species indices will be less than m_nsp - // so use volume if that is the case or area otherwise - double scalar = A; - if (cpi == spi) { - col += offset; - scalar /= A; - } else { - scalar /= m_vol; - } - // push back scaled value triplet - triplets.emplace_back(static_cast(row), static_cast(col), - scalar * it.value()); - } - } - } - // add species in this surface to the offset - offset += nk; - } + m_nv = 2 + m_nsp; // internal energy, volume, and moles of each species } void MoleReactor::getMoles(double* y) @@ -160,11 +51,6 @@ void MoleReactor::setMassFromMoles(double* y) void MoleReactor::getState(double* y) { - if (m_thermo == 0) { - throw CanteraError("MoleReactor::getState", - "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); // set the first component to the internal energy m_mass = m_thermo->density() * m_vol; y[0] = m_thermo->intEnergy_mass() * m_mass; @@ -172,9 +58,6 @@ void MoleReactor::getState(double* y) y[1] = m_vol; // set components y+2 ... y+K+2 to the moles of each species getMoles(y + m_sidx); - // set the remaining components to the surface species - // moles on walls - getSurfaceInitialConditions(y+m_nsp+m_sidx); } void MoleReactor::updateState(double* y) @@ -219,7 +102,6 @@ void MoleReactor::updateState(double* y) m_thermo->setDensity(m_mass / m_vol); } updateConnected(true); - updateSurfaceState(y + m_nsp + m_sidx); } void MoleReactor::eval(double time, double* LHS, double* RHS) @@ -227,9 +109,7 @@ void MoleReactor::eval(double time, double* LHS, double* RHS) double* dndt = RHS + m_sidx; // moles per time evalWalls(time); - m_thermo->restoreState(m_state); - - evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data()); + updateSurfaceProductionRates(); // inverse molecular weights for conversion const vector& imw = m_thermo->inverseMolecularWeights(); // volume equation @@ -290,7 +170,7 @@ size_t MoleReactor::componentIndex(const string& nm) const return 1; } try { - return speciesIndex(nm) + m_sidx; + return m_thermo->speciesIndex(nm) + m_sidx; } catch (const CanteraError&) { throw CanteraError("MoleReactor::componentIndex", "Component '{}' not found", nm); @@ -303,20 +183,7 @@ string MoleReactor::componentName(size_t k) { } else if (k == 1) { return "volume"; } else if (k >= m_sidx && k < neq()) { - k -= m_sidx; - if (k < m_thermo->nSpecies()) { - return m_thermo->speciesName(k); - } else { - k -= m_thermo->nSpecies(); - } - for (auto& S : m_surfaces) { - ThermoPhase* th = S->thermo(); - if (k < th->nSpecies()) { - return th->speciesName(k); - } else { - k -= th->nSpecies(); - } - } + return m_thermo->speciesName(k - m_sidx); } throw IndexError("MoleReactor::componentName", "component", k, m_nv); } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 58b7ee125c..54e6bf2e50 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -33,25 +33,18 @@ Reactor::Reactor(shared_ptr sol, bool clone, const string& name) m_kin = m_solution->kinetics().get(); setChemistryEnabled(m_kin->nReactions() > 0); m_vol = 1.0; // By default, the volume is set to 1.0 m^3. + m_sdot.resize(m_nsp, 0.0); + m_wdot.resize(m_nsp, 0.0); + m_nv = 3 + m_nsp; // mass, volume, internal energy, and species mass fractions } void Reactor::setDerivativeSettings(AnyMap& settings) { m_kin->setDerivativeSettings(settings); - // translate settings to surfaces - for (auto S : m_surfaces) { - S->kinetics()->setDerivativeSettings(settings); - } } void Reactor::getState(double* y) { - if (m_thermo == 0) { - throw CanteraError("Reactor::getState", - "Error: reactor is empty."); - } - m_thermo->restoreState(m_state); - // set the first component to the total mass m_mass = m_thermo->density() * m_vol; y[0] = m_mass; @@ -64,19 +57,6 @@ void Reactor::getState(double* y) // set components y+3 ... y+K+2 to the mass fractions of each species m_thermo->getMassFractions(y+3); - - // set the remaining components to the surface species - // coverages on the walls - getSurfaceInitialConditions(y + m_nsp + 3); -} - -void Reactor::getSurfaceInitialConditions(double* y) -{ - size_t loc = 0; - for (auto& S : m_surfaces) { - S->getCoverages(y + loc); - loc += S->thermo()->nSpecies(); - } } void Reactor::initialize(double t0) @@ -85,35 +65,12 @@ void Reactor::initialize(double t0) throw CanteraError("Reactor::initialize", "Reactor contents not set" " for reactor '" + m_name + "'."); } - m_thermo->restoreState(m_state); - m_sdot.resize(m_nsp, 0.0); - m_wdot.resize(m_nsp, 0.0); updateConnected(true); for (size_t n = 0; n < m_wall.size(); n++) { WallBase* W = m_wall[n]; W->initialize(); } - - m_nv = m_nsp + 3; - m_nv_surf = 0; - size_t maxnt = 0; - for (auto& S : m_surfaces) { - m_nv_surf += S->thermo()->nSpecies(); - size_t nt = S->kinetics()->nTotalSpecies(); - maxnt = std::max(maxnt, nt); - } - m_nv += m_nv_surf; - m_work.resize(maxnt); -} - -size_t Reactor::nSensParams() const -{ - size_t ns = m_sensParams.size(); - for (auto& S : m_surfaces) { - ns += S->nSensParams(); - } - return ns; } void Reactor::updateState(double* y) @@ -161,42 +118,6 @@ void Reactor::updateState(double* y) } updateConnected(true); - updateSurfaceState(y + m_nsp + 3); -} - -void Reactor::updateSurfaceState(double* y) -{ - size_t loc = 0; - for (auto& S : m_surfaces) { - S->setCoverages(y+loc); - loc += S->thermo()->nSpecies(); - } -} - -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_thermo->saveState(m_state); - - // Update the mass flow rate of connected flow devices - double time = 0.0; - if (m_net != nullptr) { - time = (timeIsIndependent()) ? m_net->time() : m_net->distance(); - } - for (size_t i = 0; i < m_outlet.size(); i++) { - m_outlet[i]->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) @@ -205,12 +126,11 @@ void Reactor::eval(double time, double* LHS, double* RHS) double* mdYdt = RHS + 3; // mass * dY/dt evalWalls(time); - m_thermo->restoreState(m_state); + updateSurfaceProductionRates(); const vector& mw = m_thermo->molecularWeights(); const double* Y = m_thermo->massFractions(); - evalSurfaces(LHS + m_nsp + 3, RHS + m_nsp + 3, m_sdot.data()); - // mass added to gas phase from surface reactions + // mass added to gas phase from surface reactions double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin()); dmdt = mdot_surf; @@ -275,35 +195,6 @@ void Reactor::evalWalls(double t) } } -void Reactor::evalSurfaces(double* LHS, double* RHS, double* sdot) -{ - fill(sdot, sdot + m_nsp, 0.0); - size_t loc = 0; // offset into ydot - - for (auto S : m_surfaces) { - Kinetics* kin = S->kinetics(); - SurfPhase* surf = S->thermo(); - - double rs0 = 1.0/surf->siteDensity(); - size_t nk = surf->nSpecies(); - double sum = 0.0; - S->restoreState(); - kin->getNetProductionRates(&m_work[0]); - for (size_t k = 1; k < nk; k++) { - RHS[loc + k] = m_work[k] * rs0 * surf->size(k); - sum -= RHS[loc + k]; - } - RHS[loc] = sum; - loc += nk; - - size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0)); - double wallarea = S->area(); - for (size_t k = 0; k < m_nsp; k++) { - sdot[k] += m_work[bulkloc + k] * wallarea; - } - } -} - vector Reactor::steadyConstraints() const { if (!energyEnabled()) { @@ -322,13 +213,7 @@ vector Reactor::steadyConstraints() const Eigen::SparseMatrix Reactor::finiteDifferenceJacobian() { - if (m_nv == 0) { - throw CanteraError("Reactor::finiteDifferenceJacobian", - "Reactor must be initialized first."); - } - // clear former jacobian elements - m_jac_trips.clear(); - + vector> trips; Eigen::ArrayXd yCurrent(m_nv); getState(yCurrent.data()); double time = (m_net != nullptr) ? m_net->time() : 0.0; @@ -359,49 +244,18 @@ Eigen::SparseMatrix Reactor::finiteDifferenceJacobian() double ydotPerturbed = rhsPerturbed[i] / lhsPerturbed[i]; double ydotCurrent = rhsCurrent[i] / lhsCurrent[i]; if (ydotCurrent != ydotPerturbed) { - m_jac_trips.emplace_back( - static_cast(i), static_cast(j), - (ydotPerturbed - ydotCurrent) / delta_y); + trips.emplace_back(static_cast(i), static_cast(j), + (ydotPerturbed - ydotCurrent) / delta_y); } } } updateState(yCurrent.data()); Eigen::SparseMatrix jac(m_nv, m_nv); - jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); + jac.setFromTriplets(trips.begin(), trips.end()); return jac; } - -void Reactor::evalSurfaces(double* RHS, double* sdot) -{ - fill(sdot, sdot + m_nsp, 0.0); - size_t loc = 0; // offset into ydot - - for (auto S : m_surfaces) { - Kinetics* kin = S->kinetics(); - SurfPhase* surf = S->thermo(); - - double rs0 = 1.0/surf->siteDensity(); - size_t nk = surf->nSpecies(); - double sum = 0.0; - S->restoreState(); - kin->getNetProductionRates(&m_work[0]); - for (size_t k = 1; k < nk; k++) { - RHS[loc + k] = m_work[k] * rs0 * surf->size(k); - sum -= RHS[loc + k]; - } - RHS[loc] = sum; - loc += nk; - - size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0)); - double wallarea = S->area(); - for (size_t k = 0; k < m_nsp; k++) { - sdot[k] += m_work[bulkloc + k] * wallarea; - } - } -} - void Reactor::addSensitivityReaction(size_t rxn) { if (!m_chem || rxn >= m_kin->nReactions()) { @@ -430,27 +284,16 @@ void Reactor::addSensitivitySpeciesEnthalpy(size_t k) SensParameterType::enthalpy}); } -size_t Reactor::speciesIndex(const string& nm) const +void Reactor::updateSurfaceProductionRates() { - // check for a gas species name - size_t k = m_thermo->speciesIndex(nm, false); - if (k != npos) { - return k; - } - - // check for a wall species - size_t offset = m_nsp; + m_sdot.assign(m_nsp, 0.0); for (auto& S : m_surfaces) { - ThermoPhase* th = S->thermo(); - k = th->speciesIndex(nm, false); - if (k != npos) { - return k + offset; - } else { - offset += th->nSpecies(); + const auto& sdot = S->surfaceProductionRates(); + size_t offset = S->kinetics()->kineticsSpeciesIndex(m_thermo->speciesName(0)); + for (size_t k = 0; k < m_nsp; k++) { + m_sdot[k] += sdot[offset + k] * S->area(); } } - throw CanteraError("Reactor::speciesIndex", - "Species '{}' not found", nm); } size_t Reactor::componentIndex(const string& nm) const @@ -465,7 +308,7 @@ size_t Reactor::componentIndex(const string& nm) const return 2; } try { - return speciesIndex(nm) + 3; + return m_thermo->speciesIndex(nm) + 3; } catch (const CanteraError&) { throw CanteraError("Reactor::componentIndex", "Component '{}' not found", nm); @@ -480,20 +323,7 @@ string Reactor::componentName(size_t k) { } else if (k == 2) { return "int_energy"; } else if (k >= 3 && k < neq()) { - k -= 3; - if (k < m_thermo->nSpecies()) { - return m_thermo->speciesName(k); - } else { - k -= m_thermo->nSpecies(); - } - for (auto& S : m_surfaces) { - ThermoPhase* th = S->thermo(); - if (k < th->nSpecies()) { - return th->speciesName(k); - } else { - k -= th->nSpecies(); - } - } + return m_thermo->speciesName(k - 3); } throw IndexError("Reactor::componentName", "component", k, m_nv); } @@ -545,9 +375,6 @@ void Reactor::applySensitivity(double* params) m_thermo->modifyOneHf298SS(p.local, p.value + params[p.global]); } } - for (auto& S : m_surfaces) { - S->setSensitivityParameters(params); - } m_thermo->invalidateCache(); if (m_kin) { m_kin->invalidateCache(); @@ -566,9 +393,6 @@ void Reactor::resetSensitivity(double* params) m_thermo->resetHf298(p.local); } } - for (auto& S : m_surfaces) { - S->resetSensitivityParameters(); - } m_thermo->invalidateCache(); if (m_kin) { m_kin->invalidateCache(); @@ -577,10 +401,6 @@ void Reactor::resetSensitivity(double* params) void Reactor::setAdvanceLimits(const double *limits) { - if (m_thermo == 0) { - throw CanteraError("Reactor::setAdvanceLimits", - "Error: reactor is empty."); - } m_advancelimits.assign(limits, limits + m_nv); // resize to zero length if no limits are set @@ -604,23 +424,6 @@ bool Reactor::getAdvanceLimits(double *limits) const void Reactor::setAdvanceLimit(const string& nm, const double limit) { size_t k = componentIndex(nm); - - if (m_thermo == 0) { - throw CanteraError("Reactor::setAdvanceLimit", - "Error: reactor is empty."); - } - if (m_nv == 0) { - if (m_net == 0) { - throw CanteraError("Reactor::setAdvanceLimit", - "Cannot set limit on a reactor that is not " - "assigned to a ReactorNet object."); - } else { - m_net->initialize(); - } - } else if (k > m_nv) { - throw CanteraError("Reactor::setAdvanceLimit", - "Index out of bounds."); - } m_advancelimits.resize(m_nv, -1.0); m_advancelimits[k] = limit; diff --git a/src/zeroD/ReactorBase.cpp b/src/zeroD/ReactorBase.cpp index 7e7af04af1..1ff243b6d2 100644 --- a/src/zeroD/ReactorBase.cpp +++ b/src/zeroD/ReactorBase.cpp @@ -7,8 +7,11 @@ #include "cantera/zeroD/FlowDevice.h" #include "cantera/zeroD/ReactorNet.h" #include "cantera/zeroD/ReactorSurface.h" +#include "cantera/zeroD/Wall.h" #include "cantera/base/Solution.h" #include "cantera/thermo/ThermoPhase.h" +#include "cantera/thermo/SurfPhase.h" +#include "cantera/kinetics/Kinetics.h" namespace Cantera { @@ -37,12 +40,10 @@ ReactorBase::ReactorBase(shared_ptr sol, bool clone, const string& nam m_solution->thermo()->addSpeciesLock(); m_thermo = m_solution->thermo().get(); m_nsp = m_thermo->nSpecies(); - m_thermo->saveState(m_state); m_enthalpy = m_thermo->enthalpy_mass(); // Needed for flow and wall interactions m_pressure = m_thermo->pressure(); // Needed for flow and wall interactions } - ReactorBase::~ReactorBase() { if (m_solution) { @@ -100,24 +101,53 @@ ReactorSurface* ReactorBase::surface(size_t n) return m_surfaces[n]; } -void ReactorBase::restoreState() { - if (!m_thermo) { - throw CanteraError("ReactorBase::restoreState", "No phase defined."); +Eigen::SparseMatrix ReactorBase::jacobian() +{ + vector> trips; + getJacobianElements(trips); + for (auto& trip : trips) { + trip = Eigen::Triplet(trip.row() - m_offset, trip.col() - m_offset, trip.value()); } - m_thermo->restoreState(m_state); + Eigen::SparseMatrix J(m_nv, m_nv); + J.setFromTriplets(trips.begin(), trips.end()); + return J; } void ReactorBase::syncState() { - m_thermo->saveState(m_state); - m_enthalpy = m_thermo->enthalpy_mass(); - m_pressure = m_thermo->pressure(); - m_mass = m_thermo->density() * m_vol; + warn_deprecated("ReactorBase::syncState", + "To be removed after Cantera 4.0. Use ReactorNet::reinitialize to indicate " + "a change in state that requires integrator reinitialization."); if (m_net) { m_net->setNeedsReinit(); } } +void ReactorBase::updateConnected(bool updatePressure) { + // save parameters needed by other connected reactors + m_enthalpy = m_thermo->enthalpy_mass(); + if (updatePressure) { + m_pressure = m_thermo->pressure(); + } + + // Update the mass flow rate of connected flow devices + double time = 0.0; + if (m_net != nullptr) { + time = (timeIsIndependent()) ? m_net->time() : m_net->distance(); + } + for (size_t i = 0; i < m_outlet.size(); i++) { + m_outlet[i]->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); + } +} + ReactorNet& ReactorBase::network() { if (m_net) { @@ -142,6 +172,26 @@ double ReactorBase::residenceTime() return mass()/mout; } +double ReactorBase::density() const +{ + return m_thermo->density(); +} + +double ReactorBase::temperature() const +{ + return m_thermo->temperature(); +} + +const double* ReactorBase::massFractions() const +{ + return m_thermo->massFractions(); +} + +double ReactorBase::massFraction(size_t k) const +{ + return m_thermo->massFraction(k); +} + FlowDevice& ReactorBase::inlet(size_t n) { return *m_inlet[n]; diff --git a/src/zeroD/ReactorFactory.cpp b/src/zeroD/ReactorFactory.cpp index bbd3030f66..ef210bd325 100644 --- a/src/zeroD/ReactorFactory.cpp +++ b/src/zeroD/ReactorFactory.cpp @@ -94,6 +94,55 @@ void ReactorFactory::deleteFactory() { s_factory = 0; } +// ---------- ReactorSurfaceFactory methods ---------- + +ReactorSurfaceFactory* ReactorSurfaceFactory::s_factory = 0; +std::mutex ReactorSurfaceFactory::s_mutex; + +ReactorSurfaceFactory* ReactorSurfaceFactory::factory() { + std::unique_lock lock(s_mutex); + if (!s_factory) { + s_factory = new ReactorSurfaceFactory; + } + return s_factory; +} + +void ReactorSurfaceFactory::deleteFactory() { + std::unique_lock lock(s_mutex); + delete s_factory; + s_factory = 0; +} + +ReactorSurfaceFactory::ReactorSurfaceFactory() +{ + reg("ReactorSurface", + [](shared_ptr surf, const vector>& reactors, + bool clone, const string& name) + { return new ReactorSurface(surf, reactors, clone, name); }); + reg("MoleReactorSurface", + [](shared_ptr surf, const vector>& reactors, + bool clone, const string& name) + { return new MoleReactorSurface(surf, reactors, clone, name); }); + reg("FlowReactorSurface", + [](shared_ptr surf, const vector>& reactors, + bool clone, const string& name) + { return new FlowReactorSurface(surf, reactors, clone, name); }); + reg("ExtensibleReactorSurface", + [](shared_ptr surf, const vector>& reactors, + bool clone, const string& name) + { return new ReactorDelegator(surf, reactors, clone, name); }); + reg("ExtensibleMoleReactorSurface", + [](shared_ptr surf, const vector>& reactors, + bool clone, const string& name) + { return new ReactorDelegator(surf, reactors, clone, name); }); + reg("ExtensibleFlowReactorSurface", + [](shared_ptr surf, const vector>& reactors, + bool clone, const string& name) + { return new ReactorDelegator(surf, reactors, clone, name); }); +} + +// ---------- free functions ---------- + shared_ptr newReactorBase( const string& model, shared_ptr phase, bool clone, const string& name) { @@ -128,7 +177,30 @@ shared_ptr newReservoir( shared_ptr newReactorSurface(shared_ptr phase, const vector>& reactors, bool clone, const string& name) { - return make_shared(phase, reactors, clone, name); + if (reactors.empty()) { + throw CanteraError("newReactorSurface", + "At least one adjacent reactor must be specified."); + } + string model = "ReactorSurface"; + if (std::dynamic_pointer_cast(reactors[0])) { + model = "FlowReactorSurface"; + } + for (const auto& r : reactors) { + if (std::dynamic_pointer_cast(r)) { + model = "MoleReactorSurface"; + break; + } + } + return shared_ptr(ReactorSurfaceFactory::factory()->create( + model, phase, reactors, clone, name)); +} + +shared_ptr newReactorSurface(const string& model, + shared_ptr phase, const vector>& reactors, + bool clone, const string& name) +{ + return shared_ptr(ReactorSurfaceFactory::factory()->create( + model, phase, reactors, clone, name)); } } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 28642a0cc7..b7fc6f4f98 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -112,27 +112,30 @@ void ReactorNet::initialize() "no reactors in network!"); } // Names of Reactors and ReactorSurfaces using each Solution; should be only one - map> solutions; + map> solutions; // Unique ReactorSurface objects. Can be attached to multiple Reactor objects - set surfaces; - m_start.assign(1, 0); - for (size_t n = 0; n < m_reactors.size(); n++) { - Reactor& r = *m_reactors[n]; + set surfaces; + m_flowDevices.clear(); + m_walls.clear(); + m_reservoirs.clear(); + m_reactors.resize(m_bulkReactors.size()); // Surfaces will be re-added + for (size_t n = 0; n < m_bulkReactors.size(); n++) { + Reactor& r = *m_bulkReactors[n]; shared_ptr bulk = r.phase(); + r.setOffset(m_nv); r.initialize(m_time); size_t nv = r.neq(); m_nv += nv; - m_start.push_back(m_nv); if (m_verbose) { writelog("Reactor {:d}: {:d} variables.\n", n, nv); writelog(" {:d} sensitivity params.\n", r.nSensParams()); } - if (r.type() == "FlowReactor" && m_reactors.size() > 1) { + if (r.type() == "FlowReactor" && m_bulkReactors.size() > 1) { throw CanteraError("ReactorNet::initialize", "FlowReactors must be used alone."); } - solutions[bulk.get()].push_back(r.name()); + solutions[bulk.get()].insert(r.name()); for (size_t i = 0; i < r.nSurfs(); i++) { if (r.surface(i)->phase()->adjacent(bulk->name()) != bulk) { throw CanteraError("ReactorNet::initialize", @@ -142,17 +145,50 @@ void ReactorNet::initialize() } surfaces.insert(r.surface(i)); } + // Discover reservoirs, flow devices, and walls connected to this reactor + for (size_t i = 0; i < r.nInlets(); i++) { + auto& inlet = r.inlet(i); + m_flowDevices.insert(&inlet); + if (inlet.in().type() == "Reservoir") { + m_reservoirs.insert(&inlet.in()); + inlet.in().setNetwork(this); + } + } + for (size_t i = 0; i < r.nOutlets(); i++) { + auto& outlet = r.outlet(i); + m_flowDevices.insert(&outlet); + if (outlet.out().type() == "Reservoir") { + m_reservoirs.insert(&outlet.out()); + outlet.out().setNetwork(this); + } + } + for (size_t i = 0; i < r.nWalls(); i++) { + auto& wall = r.wall(i); + m_walls.insert(&wall); + if (wall.left().type() == "Reservoir") { + m_reservoirs.insert(&wall.left()); + } else if (wall.right().type() == "Reservoir") { + m_reservoirs.insert(&wall.right()); + } + } } for (auto surf : surfaces) { - solutions[surf->phase().get()].push_back(surf->name()); + solutions[surf->phase().get()].insert(surf->name()); + m_reactors.push_back(surf->shared_from_this()); + surf->setOffset(m_nv); + surf->initialize(m_time); + m_nv += surf->neq(); + solutions[surf->phase().get()].insert(surf->name()); } for (auto& [soln, reactors] : solutions) { if (reactors.size() > 1) { string shared; - for (size_t i = 0; i < reactors.size() - 1; i++) { - shared += fmt::format("'{}', ", reactors[i]); + for (auto r : reactors) { + if (r != *reactors.begin()) { + shared += ", "; + } + shared += fmt::format("'{}'", r); } - shared += fmt::format("'{}'", reactors.back()); throw CanteraError("ReactorNet::initialize", "The following reactors /" " reactor surfaces are using the same Solution object: {}. Use" " independent Solution objects or set the 'clone' argument to 'true'" @@ -189,6 +225,9 @@ void ReactorNet::reinitialize() { if (m_init) { debuglog("Re-initializing reactor network.\n", m_verbose); + for (auto& reservoir : m_reservoirs) { + reservoir->updateConnected(true); + } m_integ->reinitialize(m_time, *this); if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) { checkPreconditionerSupported(); @@ -435,7 +474,7 @@ void ReactorNet::addReactor(shared_ptr reactor) reactor->type()); } - for (auto current : m_reactors) { + for (auto current : m_bulkReactors) { if (current->isOde() != r->isOde()) { throw CanteraError("ReactorNet::addReactor", "Cannot mix Reactor types using both ODEs and DAEs ({} and {})", @@ -449,7 +488,8 @@ void ReactorNet::addReactor(shared_ptr reactor) } m_timeIsIndependent = r->timeIsIndependent(); r->setNetwork(this); - m_reactors.push_back(r.get()); + m_bulkReactors.push_back(r.get()); + m_reactors.push_back(r); if (!m_integ) { m_integ.reset(newIntegrator(r->isOde() ? "CVODE" : "IDA")); // use backward differencing, with a full Jacobian computed @@ -511,19 +551,14 @@ void ReactorNet::eval(double t, double* y, double* ydot, double* p) updateState(y); m_LHS.assign(m_nv, 1); m_RHS.assign(m_nv, 0); - for (size_t n = 0; n < m_reactors.size(); n++) { - m_reactors[n]->applySensitivity(p); - m_reactors[n]->eval(t, m_LHS.data() + m_start[n], m_RHS.data() + m_start[n]); - size_t yEnd = 0; - if (n == m_reactors.size() - 1) { - yEnd = m_RHS.size(); - } else { - yEnd = m_start[n + 1]; - } - for (size_t i = m_start[n]; i < yEnd; i++) { + for (auto& R : m_reactors) { + size_t offset = R->offset(); + R->applySensitivity(p); + R->eval(t, m_LHS.data() + offset, m_RHS.data() + offset); + for (size_t i = offset; i < offset + R->neq(); i++) { ydot[i] = m_RHS[i] / m_LHS[i]; } - m_reactors[n]->resetSensitivity(p); + R->resetSensitivity(p); } checkFinite("ydot", ydot, m_nv); } @@ -532,18 +567,19 @@ void ReactorNet::evalDae(double t, double* y, double* ydot, double* p, double* r { m_time = t; updateState(y); - for (size_t n = 0; n < m_reactors.size(); n++) { - m_reactors[n]->applySensitivity(p); - m_reactors[n]->evalDae(t, y, ydot, residual); - m_reactors[n]->resetSensitivity(p); + for (auto& R : m_reactors) { + size_t offset = R->offset(); + R->applySensitivity(p); + R->evalDae(t, y + offset, ydot + offset, residual + offset); + R->resetSensitivity(p); } checkFinite("ydot", ydot, m_nv); } void ReactorNet::getConstraints(double* constraints) { - for (size_t n = 0; n < m_reactors.size(); n++) { - m_reactors[n]->getConstraints(constraints + m_start[n]); + for (auto& R : m_reactors) { + R->getConstraints(constraints + R->offset()); } } @@ -588,8 +624,8 @@ void ReactorNet::evalJacobian(double t, double* y, double* ydot, double* p, Arra void ReactorNet::updateState(double* y) { checkFinite("y", y, m_nv); - for (size_t n = 0; n < m_reactors.size(); n++) { - m_reactors[n]->updateState(y + m_start[n]); + for (auto& R : m_reactors) { + R->updateState(y + R->offset()); } } @@ -609,16 +645,16 @@ void ReactorNet::setAdvanceLimits(const double *limits) if (!m_init) { initialize(); } - for (size_t n = 0; n < m_reactors.size(); n++) { - m_reactors[n]->setAdvanceLimits(limits + m_start[n]); + for (auto& R : m_bulkReactors) { + R->setAdvanceLimits(limits + R->offset()); } } bool ReactorNet::hasAdvanceLimits() const { bool has_limit = false; - for (size_t n = 0; n < m_reactors.size(); n++) { - has_limit |= m_reactors[n]->hasAdvanceLimits(); + for (size_t n = 0; n < m_bulkReactors.size(); n++) { + has_limit |= m_bulkReactors[n]->hasAdvanceLimits(); } return has_limit; } @@ -626,23 +662,28 @@ bool ReactorNet::hasAdvanceLimits() const bool ReactorNet::getAdvanceLimits(double *limits) const { bool has_limit = false; - for (size_t n = 0; n < m_reactors.size(); n++) { - has_limit |= m_reactors[n]->getAdvanceLimits(limits + m_start[n]); + for (auto& R : m_bulkReactors) { + has_limit |= R->getAdvanceLimits(limits + R->offset()); } return has_limit; } void ReactorNet::getState(double* y) { - for (size_t n = 0; n < m_reactors.size(); n++) { - m_reactors[n]->getState(y + m_start[n]); + for (auto& R : m_reactors) { + R->getState(y + R->offset()); } } void ReactorNet::getStateDae(double* y, double* ydot) { - for (size_t n = 0; n < m_reactors.size(); n++) { - m_reactors[n]->getStateDae(y + m_start[n], ydot + m_start[n]); + // Iterate in reverse order so that surfaces will be handled first and up-to-date + // values of the surface production rates of bulk species will be available when + // bulk reactors are processed. + // TODO: Replace with view::reverse once Cantera requires C++20 + for (size_t n = m_reactors.size(); n != 0 ; n--) { + auto& R = m_reactors[n-1]; + R->getStateDae(y + R->offset(), ydot + R->offset()); } } @@ -651,7 +692,8 @@ size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor) if (!m_init) { initialize(); } - return m_start[reactor] + m_reactors[reactor]->componentIndex(component); + return m_reactors[reactor]->offset() + + m_reactors[reactor]->componentIndex(component); } string ReactorNet::componentName(size_t i) const @@ -673,7 +715,7 @@ double ReactorNet::upperBound(size_t i) const { size_t iTot = 0; size_t i0 = i; - for (auto r : m_reactors) { + for (auto r : m_bulkReactors) { if (i < r->neq()) { return r->upperBound(i); } else { @@ -688,7 +730,7 @@ double ReactorNet::lowerBound(size_t i) const { size_t iTot = 0; size_t i0 = i; - for (auto r : m_reactors) { + for (auto r : m_bulkReactors) { if (i < r->neq()) { return r->lowerBound(i); } else { @@ -700,9 +742,8 @@ double ReactorNet::lowerBound(size_t i) const } void ReactorNet::resetBadValues(double* y) { - size_t i = 0; - for (auto r : m_reactors) { - r->resetBadValues(y + m_start[i++]); + for (auto& R : m_reactors) { + R->resetBadValues(y + R->offset()); } } @@ -723,8 +764,8 @@ size_t ReactorNet::registerSensitivityParameter( void ReactorNet::setDerivativeSettings(AnyMap& settings) { // Apply given settings to all reactors - for (size_t i = 0; i < m_reactors.size(); i++) { - m_reactors[i]->setDerivativeSettings(settings); + for (auto& R : m_bulkReactors) { + R->setDerivativeSettings(settings); } } @@ -760,7 +801,8 @@ void ReactorNet::preconditionerSetup(double t, double* y, double gamma) // ensure state is up to date. updateState(y); // get the preconditioner - auto precon = m_integ->preconditioner(); + auto precon = std::dynamic_pointer_cast( + m_integ->preconditioner()); // Reset preconditioner precon->reset(); // Set gamma value for M =I - gamma*J @@ -774,15 +816,11 @@ void ReactorNet::preconditionerSetup(double t, double* y, double gamma) // update network with adjusted state updateState(yCopy.data()); // Get jacobians and give elements to preconditioners - for (size_t i = 0; i < m_reactors.size(); i++) { - Eigen::SparseMatrix rJac = m_reactors[i]->jacobian(); - for (int k=0; k::InnerIterator it(rJac, k); it; ++it) { - precon->setValue(it.row() + m_start[i], it.col() + m_start[i], - it.value()); - } - } + vector> trips; + for (auto& R : m_reactors) { + R->getJacobianElements(trips); } + precon->setFromTriplets(trips); // post reactor setup operations precon->updatePreconditioner(); } @@ -800,7 +838,7 @@ void ReactorNet::updatePreconditioner(double gamma) void ReactorNet::checkPreconditionerSupported() const { // check for non-mole-based reactors and throw an error otherwise - for (auto reactor : m_reactors) { + for (auto reactor : m_bulkReactors) { if (!reactor->preconditionerSupported()) { throw CanteraError("ReactorNet::checkPreconditionerSupported", "Preconditioning is only supported for type *MoleReactor,\n" @@ -883,6 +921,8 @@ void SteadyReactorSolver::evalJacobian(double* x0) } x0[j] = xsave; } + // Restore system to unperturbed state + m_net->updateState(x0); m_jac->updateElapsed(double(clock() - t0) / CLOCKS_PER_SEC); m_jac->incrementEvals(); diff --git a/src/zeroD/ReactorSurface.cpp b/src/zeroD/ReactorSurface.cpp index b69f77b340..92a02b6df8 100644 --- a/src/zeroD/ReactorSurface.cpp +++ b/src/zeroD/ReactorSurface.cpp @@ -48,10 +48,10 @@ ReactorSurface::ReactorSurface(shared_ptr soln, throw CanteraError("ReactorSurface::ReactorSurface", "Kinetics manager must be an InterfaceKinetics object."); } - m_kinetics = m_solution->kinetics(); + m_kinetics = std::dynamic_pointer_cast(m_solution->kinetics()); m_thermo = m_surf.get(); - m_cov.resize(m_surf->nSpecies()); - m_surf->getCoverages(m_cov.data()); + m_nsp = m_nv = m_surf->nSpecies(); + m_sdot.resize(m_kinetics->nTotalSpecies(), 0.0); } double ReactorSurface::area() const @@ -66,38 +66,85 @@ void ReactorSurface::setArea(double a) void ReactorSurface::setCoverages(const double* cov) { - copy(cov, cov + m_cov.size(), m_cov.begin()); + m_surf->setCoveragesNoNorm(cov); } void ReactorSurface::setCoverages(const Composition& cov) { m_surf->setCoveragesByName(cov); - m_surf->getCoverages(m_cov.data()); } void ReactorSurface::setCoverages(const string& cov) { m_surf->setCoveragesByName(cov); - m_surf->getCoverages(m_cov.data()); } void ReactorSurface::getCoverages(double* cov) const { - copy(m_cov.begin(), m_cov.end(), cov); + m_surf->getCoverages(cov); } -void ReactorSurface::restoreState() +void ReactorSurface::getState(double* y) { - m_surf->setTemperature(m_reactors[0]->temperature()); - m_surf->setCoveragesNoNorm(m_cov.data()); + m_surf->getCoverages(y); } -void ReactorSurface::syncState() +void ReactorSurface::initialize(double t0) { - warn_user("ReactorSurface::syncState", "Behavior changed in Cantera 3.2 for " - "consistency with ReactorBase. To set SurfPhase state from ReactorSurface " - "object, use restoreState()."); - m_surf->getCoverages(m_cov.data()); + // Sync the surface temperature and pressure to that of the first adjacent reactor + m_thermo->setState_TP(m_reactors[0]->temperature(), m_reactors[0]->pressure()); +} + +void ReactorSurface::updateState(double* y) +{ + m_surf->setCoveragesNoNorm(y); + m_thermo->setState_TP(m_reactors[0]->temperature(), m_reactors[0]->pressure()); + m_kinetics->getNetProductionRates(m_sdot.data()); +} + +void ReactorSurface::eval(double t, double* LHS, double* RHS) +{ + size_t nsp = m_surf->nSpecies(); + double rs0 = 1.0 / m_surf->siteDensity(); + double sum = 0.0; + for (size_t k = 1; k < nsp; k++) { + RHS[k] = m_sdot[k] * rs0 * m_surf->size(k); + sum -= RHS[k]; + } + RHS[0] = sum; +} + +void ReactorSurface::applySensitivity(double* params) +{ + if (!params) { + return; + } + for (auto& p : m_sensParams) { + if (p.type == SensParameterType::reaction) { + p.value = m_kinetics->multiplier(p.local); + m_kinetics->setMultiplier(p.local, p.value*params[p.global]); + } else if (p.type == SensParameterType::enthalpy) { + m_thermo->modifyOneHf298SS(p.local, p.value + params[p.global]); + } + } + m_thermo->invalidateCache(); + m_kinetics->invalidateCache(); +} + +void ReactorSurface::resetSensitivity(double* params) +{ + if (!params) { + return; + } + for (auto& p : m_sensParams) { + if (p.type == SensParameterType::reaction) { + m_kinetics->setMultiplier(p.local, p.value); + } else if (p.type == SensParameterType::enthalpy) { + m_thermo->resetHf298(p.local); + } + } + m_thermo->invalidateCache(); + m_kinetics->invalidateCache(); } void ReactorSurface::addSensitivityReaction(size_t rxn) @@ -112,19 +159,146 @@ void ReactorSurface::addSensitivityReaction(size_t rxn) SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction}); } -void ReactorSurface::setSensitivityParameters(const double* params) +size_t ReactorSurface::componentIndex(const string& nm) const { - for (auto& p : m_sensParams) { - p.value = m_kinetics->multiplier(p.local); - m_kinetics->setMultiplier(p.local, p.value*params[p.global]); + return m_surf->speciesIndex(nm); +} + +string ReactorSurface::componentName(size_t k) +{ + return m_surf->speciesName(k); +} + +// ------ MoleReactorSurface methods ------ + +void MoleReactorSurface::initialize(double t0) +{ + ReactorSurface::initialize(t0); + m_cov_tmp.resize(m_nsp); + m_f_energy.resize(m_kinetics->nTotalSpecies(), 0.0); + m_f_species.resize(m_kinetics->nTotalSpecies(), 0.0); + m_kin2net.resize(m_kinetics->nTotalSpecies(), -1); + m_kin2reactor.resize(m_kinetics->nTotalSpecies(), nullptr); + + for (size_t k = 0; k < m_nsp; k++) { + m_kin2net[k] = static_cast(m_offset + k); + } + for (auto R : m_reactors) { + size_t nsp = R->phase()->thermo()->nSpecies(); + size_t k0 = m_kinetics->speciesOffset(*R->phase()->thermo()); + int offset = static_cast(R->offset() + R->speciesOffset()); + for (size_t k = 0; k < nsp; k++) { + m_kin2net[k + k0] = offset + k; + m_kin2reactor[k + k0] = R; + } } } -void ReactorSurface::resetSensitivityParameters() +void MoleReactorSurface::getState(double* y) { - for (auto& p : m_sensParams) { - m_kinetics->setMultiplier(p.local, p.value); + m_surf->getCoverages(y); + double totalSites = m_surf->siteDensity() * m_area; + for (size_t k = 0; k < m_nsp; k++) { + y[k] *= totalSites / m_surf->size(k); } } +void MoleReactorSurface::updateState(double* y) +{ + std::copy(y, y + m_nsp, m_cov_tmp.data()); + double totalSites = m_surf->siteDensity() * m_area; + for (size_t k = 0; k < m_nsp; k++) { + m_cov_tmp[k] *= m_surf->size(k) / totalSites; + } + m_surf->setCoveragesNoNorm(m_cov_tmp.data()); + m_thermo->setState_TP(m_reactors[0]->temperature(), m_reactors[0]->pressure()); + m_kinetics->getNetProductionRates(m_sdot.data()); +} + +void MoleReactorSurface::eval(double t, double* LHS, double* RHS) +{ + for (size_t k = 0; k < m_nsp; k++) { + RHS[k] = m_sdot[k] * m_area / m_surf->size(k); + } +} + +void MoleReactorSurface::getJacobianElements(vector>& trips) +{ + auto sdot_ddC = m_kinetics->netProductionRates_ddCi(); + for (auto R : m_reactors) { + double f_species; + size_t nsp_R = R->phase()->thermo()->nSpecies(); + size_t k0 = m_kinetics->speciesOffset(*R->phase()->thermo()); + R->getJacobianScalingFactors(f_species, &m_f_energy[k0]); + std::fill(&m_f_species[k0], &m_f_species[k0 + nsp_R], m_area * f_species); + } + std::fill(&m_f_species[0], &m_f_species[m_nsp], 1.0); // surface species + for (int k = 0; k < sdot_ddC.outerSize(); k++) { + int col = m_kin2net[k]; + if (col == -1) { + continue; + } + for (Eigen::SparseMatrix::InnerIterator it(sdot_ddC, k); it; ++it) { + int row = m_kin2net[it.row()]; + if (row == -1 || it.value() == 0.0) { + continue; + } + ReactorBase* R = m_kin2reactor[it.row()]; + trips.emplace_back(row, col, it.value() * m_f_species[k]); + if (m_f_energy[it.row()] != 0.0) { + trips.emplace_back(R->offset(), col, + it.value() * m_f_energy[it.row()] * m_f_species[k]); + } + } + } +} + +// ------ FlowReactorSurface methods ------ + +FlowReactorSurface::FlowReactorSurface(shared_ptr soln, + const vector>& reactors, + bool clone, + const string& name) + : ReactorSurface(soln, reactors, clone, name) +{ + m_area = -1.0; // default to perimeter of cylindrical reactor +} + +double FlowReactorSurface::area() const { + if (m_area > 0) { + return m_area; + } + + // Assuming a cylindrical cross section, P = 2 * pi * r and A = pi * r^2, so + // P(A) = 2 * sqrt(pi * A) + return 2.0 * sqrt(Pi * m_reactors[0]->area()); +} + +void FlowReactorSurface::evalDae(double t, double* y, double* ydot, double* residual) +{ + size_t nsp = m_surf->nSpecies(); + double sum = y[0]; + for (size_t k = 1; k < nsp; k++) { + residual[k] = m_sdot[k]; + sum += y[k]; + } + residual[0] = sum - 1.0; +} + +void FlowReactorSurface::getStateDae(double* y, double* ydot) +{ + // Advance the surface to steady state to get consistent initial coverages + m_kinetics->advanceCoverages(100.0, m_ss_rtol, m_ss_atol, 0, m_max_ss_steps, + m_max_ss_error_fails); + getCoverages(y); + // Update the values in m_sdot that are needed by the adjacent FlowReactor + updateState(y); +} + +void FlowReactorSurface::getConstraints(double* constraints) +{ + // The species coverages are algebraic constraints + std::fill(constraints, constraints + m_nsp, 1.0); +} + } diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 736e95ee8a..f1627f6d46 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -933,7 +933,7 @@ def test_set_initial_time(self): assert p1a == approx(p1b) assert p2a == approx(p2b) - def test_reinitialize(self): + def test_reinitialize(self, allow_deprecated): self.make_reactors(T1=300, T2=1000) self.add_wall(U=200, A=1.0) self.net.advance(1.0) @@ -941,6 +941,7 @@ def test_reinitialize(self): T2a = self.r2.T self.r1.phase.TD = 300, None + # Deprecated; After Cantera 4.0, replace with net.reinitialize() self.r1.syncState() self.r2.phase.TD = 1000, None @@ -955,6 +956,13 @@ def test_reinitialize(self): assert T1a == approx(T1b) assert T2a == approx(T2b) + # TODO: Remove after Cantera 4.0 when syncState is removed + def test_syncState_deprecated(self): + self.make_reactors(n_reactors=1) + self.net.advance(0.1) + with pytest.raises(ct.CanteraError): + self.r1.syncState() + def test_reservoir_sync(self): self.make_reactors(T1=800, n_reactors=1) self.gas1.TP = 900, ct.one_atm @@ -963,7 +971,6 @@ def test_reservoir_sync(self): self.net.advance(1.0) assert self.r1.T == approx(872.099, rel=1e-3) reservoir.phase.TP = 700, ct.one_atm - reservoir.syncState() self.net.reinitialize() self.net.advance(2.0) assert self.r1.T == approx(747.27, rel=1e-3) @@ -1452,16 +1459,15 @@ def test_reactor_surface_type(self): def test_component_index(self): self.create_reactors(add_surf=True) - for (gas,net,iface,r) in ((self.gas1, self.net1, self.interface1, self.r1), - (self.gas2, self.net2, self.interface2, self.r2)): + for (gas,net,iface,r,s) in ((self.gas1, self.net1, self.interface1, self.r1, self.surf1), + (self.gas2, self.net2, self.interface2, self.r2, self.surf2)): net.step() N0 = net.n_vars - gas.n_species - iface.n_species - N1 = net.n_vars - iface.n_species for i, name in enumerate(gas.species_names): assert i + N0 == r.component_index(name) for i, name in enumerate(iface.species_names): - assert i + N1 == r.component_index(name) + assert i == s.component_index(name) def test_component_names(self): self.create_reactors(add_surf=True) @@ -1818,11 +1824,11 @@ def test_catalytic_surface(self): porosity = 0.3 velocity = 0.4 / 60 mdot = velocity * gas.density * r.area * porosity - r.surface_area_to_volume_ratio = porosity * 1e5 r.mass_flow_rate = mdot r.energy_enabled = False - rsurf = ct.ReactorSurface(surf, r) + rsurf = ct.FlowReactorSurface(surf, r) + rsurf.area = 1e5 * porosity * r.area sim = ct.ReactorNet([r]) kCH4 = gas.species_index('CH4') @@ -1856,12 +1862,19 @@ def test_component_names(self): sim = ct.ReactorNet([r]) sim.initialize() - assert r.n_vars == 4 + gas.n_species + surf.n_species - assert sim.n_vars == r.n_vars + assert r.n_vars == 4 + gas.n_species + assert rsurf.n_vars == surf.n_species + assert sim.n_vars == r.n_vars + rsurf.n_vars for i in range(r.n_vars): name = r.component_name(i) assert r.component_index(name) == i + assert name in sim.component_name(i) + + for i in range(rsurf.n_vars): + name = rsurf.component_name(i) + assert rsurf.component_index(name) == i + assert name in sim.component_name(i + r.n_vars) with pytest.raises(ct.CanteraError, match="Component 'spam' not found"): r.component_index('spam') @@ -1869,6 +1882,7 @@ def test_component_names(self): with pytest.raises(ct.CanteraError, match="outside valid range"): r.component_name(200) + class TestFlowReactor2: def import_phases(self): surf = ct.Interface('SiF4_NH3_mec.yaml', 'SI3N4') @@ -1877,9 +1891,9 @@ def import_phases(self): def make_reactors(self, gas, surf): r = ct.FlowReactor(gas) r.area = 1e-4 - r.surface_area_to_volume_ratio = 5000 r.mass_flow_rate = 0.02 - rsurf = ct.ReactorSurface(surf, r) + rsurf = ct.FlowReactorSurface(surf, r) + rsurf.area = 5000 * r.area sim = ct.ReactorNet([r]) return r, rsurf, sim @@ -2082,11 +2096,10 @@ def test_reinitialization(self): # Reset the reactor to the same initial state r.phase.TPX = 1700, 4000, 'NH3:1.0, SiF4:0.4' - surf.TP = 1700, 4000 + rsurf.phase.TP = 1700, 4000 r.mass_flow_rate = 0.01 - r.syncState() - sim.initial_time = 0. + sim.reinitialize() sim.advance(0.6) Ygas2 = r.phase.Y cov2 = rsurf.phase.coverages @@ -2102,26 +2115,25 @@ def test_initial_condition_tolerances(self): r, rsurf, sim = self.make_reactors(gas, surf) # With tight tolerances, some error test failures should be expected - r.inlet_surface_atol = 1e-14 - r.inlet_surface_rtol = 1e-20 - r.inlet_surface_max_error_failures = 1 + rsurf.initial_atol = 1e-14 + rsurf.initial_rtol = 1e-20 + rsurf.initial_max_error_failures = 1 with pytest.raises(ct.CanteraError, match="error test failed repeatedly"): sim.initialize() # With few steps allowed, won't be able to reach steady-state - r.inlet_surface_max_error_failures = 10 - r.inlet_surface_max_steps = 200 + rsurf.initial_max_error_failures = 10 + rsurf.initial_max_steps = 200 with pytest.raises(ct.CanteraError, match="Maximum number of timesteps"): sim.initialize() # Relaxing the tolerances will allow the integrator to reach steady-state # in fewer timesteps surf.coverages = np.ones(surf.n_species) - r.inlet_surface_atol = 0.001 - r.inlet_surface_rtol = 0.001 + rsurf.initial_atol = 0.001 + rsurf.initial_rtol = 0.001 sim.initialize() - def test_Si3N4_deposition_regression(): # Regression test based on silicon nitride deposition example given in # 1D_pfr_surfchem.py with values published in Sandia Report SAND-96-8211 @@ -2384,7 +2396,7 @@ def test_sensitivities2(self): Ns = r1.component_index(gas1.species_name(0)) # Index of first variable corresponding to r2 - K2 = Ns + gas1.n_species + interface.n_species + K2 = Ns + gas1.n_species # Constant volume should generate zero sensitivity coefficient assert S[1,:] == approx(np.zeros(2)) @@ -3009,8 +3021,6 @@ def __init__(self, *args, neighbor, **kwargs): self.v_wall = 0 self.k_wall = 1e-5 self.neighbor = neighbor - - def after_initialize(self, t0): self.n_vars += 1 self.i_wall = self.n_vars - 1 @@ -3105,6 +3115,13 @@ def replace_component_index(self, name): r2.component_index("H2") assert r2.component_index("succeed") == 0 + class DummyReactor3(ct.ExtensibleReactor): + def before_foobar(self, t): + pass + + with pytest.raises(ValueError, match="'foobar' is not a delegatable method"): + DummyReactor3(self.gas) + def test_delegate_throws(self): class TestException(Exception): pass @@ -3135,23 +3152,19 @@ def test_misc(self): class DummyReactor(ct.ExtensibleReactor): def __init__(self, gas, *args, **kwargs): super().__init__(gas, *args, **kwargs) - self.sync_calls = 0 - def after_species_index(self, name): + def after_component_index(self, name): # This will cause returned species indices to be higher by 5 than they # would be otherwise - return 5 - - def before_sync_state(self): - self.sync_calls += 1 + if name in self.phase.species_names: + return 5 + else: + return 0 r = DummyReactor(self.gas) net = ct.ReactorNet([r]) assert r.component_index("H2") == 5 + 3 + self.gas.species_index("H2") - r.syncState() - net.advance(1) - r.syncState() - assert r.sync_calls == 2 + assert r.component_index("int_energy") == 2 def test_RHS_LHS(self): # set initial state @@ -3241,41 +3254,40 @@ def test_with_surface(self): kHs = surf.species_index("H(S)") kPts = surf.species_index("PT(S)") kH2 = gas.species_index("H2") + kH2_kin = surf.kinetics_species_index("H2") kO2 = gas.species_index("O2") - class SurfReactor(ct.ExtensibleIdealGasReactor): - def replace_eval_surfaces(self, LHS, RHS, sdot): - site_density = self.surfaces[0].phase.site_density - sdot[:] = 0.0 - LHS[:] = 1.0 - RHS[:] = 0.0 - C = self.phase.concentrations - theta = self.surfaces[0].coverages - - # Replace actual reactions with simple absorption of H2 -> H(S) - rop = 1e-4 * C[kH2] * theta[kPts] - RHS[kHs] = 2 * rop / site_density - RHS[kPts] = - 2 * rop / site_density - sdot[kH2] = - rop * self.surfaces[0].area - - def replace_get_surface_initial_conditions(self, y): + class CustomSurface(ct.ExtensibleReactorSurface): + def replace_get_state(self, y): y[:] = 0 y[kPts] = 1 - def replace_update_surface_state(self, y): + def replace_update_state(self, y): # this is the same thing the original method does - self.surfaces[0].coverages = y + self.phase.set_unnormalized_coverages(y) + + # Replace actual reactions with simple absorption of H2 -> H(S) + C = self.reactors[0].phase.concentrations + theta = self.phase.coverages + self.rop = 1e-4 * C[kH2] * theta[kPts] + self.surface_production_rates[:] = 0.0 + self.surface_production_rates[kH2_kin] = - self.rop - r1 = SurfReactor(gas) + def replace_eval(self, t, LHS, RHS): + site_density = self.phase.site_density + RHS[kHs] = 2 * self.rop / site_density + RHS[kPts] = - 2 * self.rop / site_density + + r1 = ct.IdealGasReactor(gas) r1.volume = 1e-6 # 1 cm^3 r1.energy_enabled = False - rsurf = ct.ReactorSurface(surf, r=r1, A=0.01) + rsurf = CustomSurface(surf, r=r1, A=0.01) net = ct.ReactorNet([r1]) Hweight = ct.Element("H").weight total_sites = rsurf.area * surf.site_density def masses(): mass_H = (r1.phase.elemental_mass_fraction("H") * r1.mass + - total_sites * r1.surfaces[0].phase["H(s)"].X * Hweight) + total_sites * rsurf.phase["H(S)"].X * Hweight) mass_O = r1.phase.elemental_mass_fraction("O") * r1.mass return mass_H, mass_O @@ -3292,8 +3304,8 @@ def masses(): assert r1.phase.P == approx(647.56016304) assert r1.phase.X[kH2] == approx(0.4784268406) assert r1.phase.X[kO2] == approx(0.5215731594) - assert r1.surfaces[0].phase.X[kHs] == approx(0.3665198138) - assert r1.surfaces[0].phase.X[kPts] == approx(0.6334801862) + assert rsurf.phase.X[kHs] == approx(0.3665198138) + assert rsurf.phase.X[kPts] == approx(0.6334801862) def test_interactions(self): # Reactors connected by a movable, H2-permeable surface