From 205567a99b1af4964af060ab6d28f22091fe12e9 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Tue, 17 Oct 2023 10:33:12 -0700 Subject: [PATCH] Build jacobian structure defined for walls and flow devices --- include/cantera/zeroD/FlowDevice.h | 41 ++++++++++++ .../zeroD/IdealGasConstPressureMoleReactor.h | 12 ++++ include/cantera/zeroD/IdealGasMoleReactor.h | 12 ++++ include/cantera/zeroD/MoleReactor.h | 6 ++ include/cantera/zeroD/Reactor.h | 26 +++++++- include/cantera/zeroD/ReactorBase.h | 58 +++++++++++++++++ include/cantera/zeroD/ReactorNet.h | 18 ++++++ include/cantera/zeroD/Wall.h | 49 ++++++++++++++ interfaces/cython/cantera/reactor.pxd | 2 + interfaces/cython/cantera/reactor.pyx | 26 ++++++++ .../IdealGasConstPressureMoleReactor.cpp | 26 ++++++++ src/zeroD/IdealGasMoleReactor.cpp | 26 ++++++++ src/zeroD/Reactor.cpp | 33 ++++++++++ src/zeroD/ReactorNet.cpp | 55 ++++++++++++++++ src/zeroD/Wall.cpp | 64 +++++++++++++++++++ test/python/test_reactor.py | 64 ++++++++++++++++--- 16 files changed, 508 insertions(+), 10 deletions(-) diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 463cdb7a9c0..9590d374988 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -9,6 +9,7 @@ #include "cantera/base/ct_defs.h" #include "cantera/base/global.h" #include "cantera/base/ctexceptions.h" +#include "cantera/numerics/eigen_sparse.h" namespace Cantera { @@ -114,7 +115,47 @@ class FlowDevice m_time = time; } + /*! Build the Jacobian terms specific to the flow device for the given connected reactor. + * @param r a pointer to the calling reactor + * @param jacVector a vector of triplets to be added to the jacobian for the reactor + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { + throw NotImplementedError(type() + "::buildReactorJacobian"); + } + + /*! Build the Jacobian terms specific to the flow device for the network. These terms + * will be adjusted to the networks indexing system outside of the reactor. + * @param jacVector a vector of triplets to be added to the jacobian for the reactor + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual void buildNetworkJacobian(vector>& jacVector) { + throw NotImplementedError(type() + "::buildNetworkJacobian"); + } + + /*! Specify the jacobian terms have been calculated and should not be recalculated. + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + void jacobianCalculated() { m_jac_calculated = true; }; + + /*! Specify that jacobian terms have not been calculated and should be recalculated. + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + void jacobianNotCalculated() { m_jac_calculated = false; }; + protected: + //! a variable to switch on and off so calculations are not doubled by the calling + //! reactor or network + bool m_jac_calculated = false; + double m_mdot = Undef; //! Function set by setPressureFunction; used by updateMassFlowRate diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 3c6dc207213..3064a641134 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -49,6 +49,18 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor bool preconditionerSupported() const override { return true; }; + double temperatureDerivative() override { return 1; }; + + double temperatureRadiationDerivative() override { + return 4 * std::pow(temperature(), 3); + } + + double moleDerivative(size_t index) override; + + double moleRadiationDerivative(size_t index) override; + + size_t speciesOffset() const override { return m_sidx; }; + protected: vector m_hk; //!< Species molar enthalpies }; diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index 48570556876..82cc9f70fbc 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -45,6 +45,18 @@ class IdealGasMoleReactor : public MoleReactor bool preconditionerSupported() const override {return true;}; + double temperatureDerivative() override { return 1; }; + + double temperatureRadiationDerivative() override { + return 4 * std::pow(temperature(), 3); + } + + double moleDerivative(size_t index) override; + + double moleRadiationDerivative(size_t index) override; + + size_t speciesOffset() const override { return m_sidx; }; + protected: vector m_uk; //!< Species molar internal energies }; diff --git a/include/cantera/zeroD/MoleReactor.h b/include/cantera/zeroD/MoleReactor.h index f543d7d7838..cddf099f388 100644 --- a/include/cantera/zeroD/MoleReactor.h +++ b/include/cantera/zeroD/MoleReactor.h @@ -7,6 +7,7 @@ #define CT_MOLEREACTOR_H #include "Reactor.h" +#include "cantera/zeroD/Wall.h" namespace Cantera { @@ -38,6 +39,8 @@ class MoleReactor : public Reactor string componentName(size_t k) override; + size_t energyIndex() const override { return m_eidx; }; + protected: //! For each surface in the reactor, update vector of triplets with all relevant //! surface jacobian derivatives of species with respect to species @@ -60,6 +63,9 @@ class MoleReactor : public Reactor //! const value for the species start index const size_t m_sidx = 2; + + //! index of state variable associated with energy + const size_t m_eidx = 0; }; } diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 67229c006f8..b6af3c2ee2a 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -232,7 +232,27 @@ class Reactor : public ReactorBase throw NotImplementedError(type() + "::buildJacobian"); } - // virtual void jacobian() + //! Calculate the Jacobian of a Reactor specialization for wall contributions. + //! @param jacVector vector where jacobian triplets are added + //! @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 void buildWallJacobian(vector>& jacVector); + + //! Calculate flow contributions to the Jacobian of a Reactor specialization. + //! @param jac_vector vector where jacobian triplets are added + //! @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 void buildFlowJacobian(vector>& jacVector); //! Calculate the reactor-specific Jacobian using a finite difference method. //! @@ -325,6 +345,10 @@ class Reactor : public ReactorBase //! Vector of triplets representing the jacobian vector> m_jac_trips; + //! Boolean to skip walls in jacobian + bool m_jac_skip_walls = false; + //! Boolean to skip flow devices in jacobian + bool m_jac_skip_flow_devices = false; }; } diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index f5f3a35fe00..43b0165e42a 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -248,10 +248,68 @@ class ReactorBase //! Set the ReactorNet that this reactor belongs to. void setNetwork(ReactorNet* net); + /*! Calculate the derivative of temperature with respect to the temperature in the + * heat transfer equation based on the reactor specific equation of state. + * This function should also transform the state of the derivative to that + * appropriate for the jacobian's state/ + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual double temperatureDerivative() { + throw NotImplementedError("Reactor::temperatureDerivative"); + } + + /*! Calculate the derivative of temperature with respect to the temperature in the + * heat transfer radiation equation based on the reactor specific equation of state. + * This function should also transform the state of the derivative to that + * appropriate for the jacobian's state/ + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual double temperatureRadiationDerivative() { + throw NotImplementedError("Reactor::temperatureRadiationDerivative"); + } + + /*! Calculate the derivative of T with respect to the ith species in the heat + * transfer equation based on the reactor specific equation of state. + * @param index index of the species the derivative is with respect too + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual double moleDerivative(size_t index) { + throw NotImplementedError("Reactor::moleDerivative"); + } + + /*! Calculate the derivative of T with respect to the ith species in the heat + * transfer radiation equation based on the reactor specific equation of state. + * @param index index of the species the derivative is with respect too + * @warning This function is an experimental part of the %Cantera API and may be changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual double moleRadiationDerivative(size_t index) { + throw NotImplementedError("Reactor::moleRadiationDerivative"); + } + + //! Return the index associated with energy of the system + virtual size_t energyIndex() const { return m_eidx; }; + + //! Return the offset between species and state variables + virtual size_t speciesOffset() const { return m_sidx; }; + protected: //! Number of homogeneous species in the mixture size_t m_nsp = 0; + //! species offset in the state vector + const size_t m_sidx = 3; + + //! index of state variable associated with energy + const size_t m_eidx = 1; + ThermoPhase* m_thermo = nullptr; double m_vol = 1.0; //!< Current volume of the reactor [m^3] double m_enthalpy = 0.0; //!< Current specific enthalpy of the reactor [J/kg] diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 55a48b43531..14367184805 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -236,6 +236,20 @@ class ReactorNet : public FuncEval //! reactor network. size_t globalComponentIndex(const string& component, size_t reactor=0); + //! Return the index corresponding to the start of the reactor specific state + //! vector in the reactor with index *reactor* in the global state vector for the + //! reactor network. + size_t globalStartIndex(Reactor* curr_reactor) { + for (size_t i = 0; i < m_reactors.size(); i++) { + if (curr_reactor == m_reactors[i]) { + return m_start[i]; + } + } + throw CanteraError("ReactorNet::globalStartIndex: ", + curr_reactor->name(), " not found in network."); + return npos; + } + //! Return the name of the i-th component of the global state vector. The //! name returned includes both the name of the reactor and the specific //! component, for example `'reactor1: CH4'`. @@ -396,6 +410,10 @@ class ReactorNet : public FuncEval //! "left hand side" of each governing equation vector m_LHS; vector m_RHS; + + //! derivative settings + bool m_jac_skip_walls = false; + bool m_jac_skip_flow_devices = false; }; } diff --git a/include/cantera/zeroD/Wall.h b/include/cantera/zeroD/Wall.h index e3e1ec8cb5d..d240ff4f7e6 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -8,6 +8,7 @@ #include "cantera/base/ctexceptions.h" #include "cantera/zeroD/ReactorBase.h" +#include "cantera/numerics/eigen_sparse.h" namespace Cantera { @@ -115,6 +116,46 @@ class WallBase m_time = time; } + /*! Build the Jacobian terms specific to the flow device for the given connected reactor. + * @param r a pointer to the calling reactor + * @param jacVector a vector of triplets to be added to the jacobian for the reactor + * @warning This function is an experimental part of the %Cantera API and may be + * changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { + throw NotImplementedError("WallBase::buildReactorJacobian"); + } + + /*! Build the Jacobian terms specific to the flow device for the network. These terms + * will be adjusted to the networks indexing system outside of the reactor. + * @param jacVector a vector of triplets to be added to the jacobian for the reactor + * @warning This function is an experimental part of the %Cantera API and may be + * changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + virtual void buildNetworkJacobian(vector>& jacVector) { + throw NotImplementedError("WallBase::buildNetworkJacobian"); + } + + /*! Specify the jacobian terms have been calculated and should not be recalculated. + * @warning This function is an experimental part of the %Cantera API and may be + * changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + void jacobianCalculated() { m_jac_calculated = true; }; + + /*! Specify that jacobian terms have not been calculated and should be recalculated. + * @warning This function is an experimental part of the %Cantera API and may be + * changed + * or removed without notice. + * @since New in %Cantera 3.0. + */ + void jacobianNotCalculated() { m_jac_calculated = false; }; + protected: ReactorBase* m_left = nullptr; ReactorBase* m_right = nullptr; @@ -123,6 +164,10 @@ class WallBase double m_time = 0.0; double m_area = 1.0; + + //! a variable to switch on and off so calculations are not doubled by the calling + //! reactor or network + bool m_jac_calculated = false; }; //! Represents a wall between between two ReactorBase objects. @@ -255,6 +300,10 @@ class Wall : public WallBase return m_k; } + virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) override; + + virtual void buildNetworkJacobian(vector>& jacVector) override; + protected: //! expansion rate coefficient diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index bedc6e7d1d9..a528b7b6e72 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -202,6 +202,8 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": void setPreconditioner(shared_ptr[CxxPreconditionerBase] preconditioner) void setDerivativeSettings(CxxAnyMap&) CxxAnyMap solverStats() except +translate_exception + CxxSparseMatrix jacobian() except +translate_exception + CxxSparseMatrix finiteDifferenceJacobian() except +translate_exception cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera": cdef cppclass CxxReactorAccessor "Cantera::ReactorAccessor": diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 1cded0cb5a1..d0aa5270864 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -1797,3 +1797,29 @@ cdef class ReactorNet: """ def __set__(self, settings): self.net.setDerivativeSettings(py_to_anymap(settings)) + + property jacobian: + """ + Get the system Jacobian or an approximation thereof. + + **Warning**: Depending on the particular implementation, this may return an + approximate Jacobian intended only for use in forming a preconditioner for + iterative solvers, excluding terms that would generate a fully-dense Jacobian. + + **Warning**: This method is an experimental part of the Cantera API and may be + changed or removed without notice. + """ + def __get__(self): + return get_from_sparse(self.net.jacobian(), self.n_vars, self.n_vars) + + property finite_difference_jacobian: + """ + Get the system Jacobian, calculated using a finite difference method. + + **Warning:** this property is an experimental part of the Cantera API and + may be changed or removed without notice. + """ + def __get__(self): + return get_from_sparse(self.net.finiteDifferenceJacobian(), + self.n_vars, self.n_vars) + diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index 912a7736171..59c807464c5 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -203,6 +203,7 @@ void IdealGasConstPressureMoleReactor::buildJacobian( jacVector.emplace_back(static_cast(j), 0, (ydotPerturbed - ydotCurrent) / deltaTemp); } + // d T_dot/dnj // allocate vectors for whole system Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize); @@ -228,7 +229,13 @@ void IdealGasConstPressureMoleReactor::buildJacobian( jacVector.emplace_back(0, static_cast(j + m_sidx), (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom); } + + // build wall jacobian + buildWallJacobian(jacVector); } + + // build flow jacobian + buildFlowJacobian(jacVector); } size_t IdealGasConstPressureMoleReactor::componentIndex(const string& nm) const @@ -266,4 +273,23 @@ string IdealGasConstPressureMoleReactor::componentName(size_t k) { "Index is out of bounds."); } +double IdealGasConstPressureMoleReactor::moleDerivative(size_t index) +{ + // derivative of temperature transformed by ideal gas law + vector moles(m_nsp); + getMoles(moles.data()); + double dTdni = pressure() * m_vol / GasConstant / std::accumulate(moles.begin(), moles.end(), 0.0); + return dTdni; +} + +double IdealGasConstPressureMoleReactor::moleRadiationDerivative(size_t index) +{ + // derivative of temperature transformed by ideal gas law + vector moles(m_nsp); + getMoles(moles.data()); + double dT4dni = std::pow(pressure() * m_vol / GasConstant, 4); + dT4dni *= std::pow(1 / std::accumulate(moles.begin(), moles.end(), 0.0), 5); + return dT4dni; +} + } diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index fd51a1c460b..9af258b4719 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -214,6 +214,7 @@ void IdealGasMoleReactor::buildJacobian(vector>& jacVecto jacVector.emplace_back(static_cast(j), 0, (ydotPerturbed - ydotCurrent) / deltaTemp); } + // d T_dot/dnj Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize); Eigen::VectorXd internal_energy = Eigen::VectorXd::Zero(ssize); @@ -243,7 +244,32 @@ void IdealGasMoleReactor::buildJacobian(vector>& jacVecto jacVector.emplace_back(0, static_cast(j + m_sidx), (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom); } + + // build wall jacobian + buildWallJacobian(jacVector); } + + // build flow jacobian + buildFlowJacobian(jacVector); +} + +double IdealGasMoleReactor::moleDerivative(size_t index) +{ + // derivative of temperature transformed by ideal gas law + vector moles(m_nsp); + getMoles(moles.data()); + double dTdni = pressure() * m_vol / GasConstant / std::accumulate(moles.begin(), moles.end(), 0.0); + return dTdni; +} + +double IdealGasMoleReactor::moleRadiationDerivative(size_t index) +{ + // derivative of temperature transformed by ideal gas law + vector moles(m_nsp); + getMoles(moles.data()); + double dT4dni = std::pow(pressure() * m_vol / GasConstant, 4); + dT4dni *= std::pow(1 / std::accumulate(moles.begin(), moles.end(), 0.0), 5); + return dT4dni; } } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 49ae58844d1..5dfd31433b9 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -34,6 +34,17 @@ void Reactor::setDerivativeSettings(AnyMap& settings) for (auto S : m_surfaces) { S->kinetics()->setDerivativeSettings(settings); } + + // set reactor settings + bool force = settings.empty(); + if (force || settings.hasKey("skip-walls")) { + m_jac_skip_walls = settings.getBool("skip-walls", + false); + } + if (force || settings.hasKey("skip-flow-devices")) { + m_jac_skip_flow_devices = settings.getBool("skip-flow-devices", + false); + } } void Reactor::setKineticsMgr(Kinetics& kin) @@ -589,4 +600,26 @@ void Reactor::setAdvanceLimit(const string& nm, const double limit) } } +void Reactor:: buildWallJacobian(vector>& jacVector) +{ + if (!m_jac_skip_walls) { + for (size_t i = 0; i < m_wall.size(); i++) { + m_wall[i]->buildReactorJacobian(this, jacVector); + } + } +} + +void Reactor::buildFlowJacobian(vector>& jacVector) +{ + if (!m_jac_skip_flow_devices) { + for (size_t i = 0; i < m_outlet.size(); i++) { + m_outlet[i]->buildReactorJacobian(this, jacVector); + } + + for (size_t i = 0; i buildReactorJacobian(this, jacVector); + } + } +} + } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 64c2c742106..4261bfe03c9 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -496,6 +496,16 @@ void ReactorNet::setDerivativeSettings(AnyMap& settings) for (size_t i = 0; i < m_reactors.size(); i++) { m_reactors[i]->setDerivativeSettings(settings); } + // set network settings + bool force = settings.empty(); + if (force || settings.hasKey("skip-walls")) { + m_jac_skip_walls = settings.getBool("skip-walls", + false); + } + if (force || settings.hasKey("skip-flow-devices")) { + m_jac_skip_flow_devices = settings.getBool("skip-flow-devices", + false); + } } AnyMap ReactorNet::solverStats() const @@ -583,6 +593,26 @@ void ReactorNet::buildJacobian(vector>& jacVector) if (! m_init) { initialize(); } + // loop through and set connectors not found + for (auto r : m_reactors) { + // walls + if (!m_jac_skip_walls) { + for (size_t i = 0; i < r->nWalls(); i++) { + r->wall(i).jacobianNotCalculated(); + } + } + // flow devices + if (!m_jac_skip_flow_devices) { + // outlets + for (size_t i = 0; i < r->nOutlets(); i++) { + r->outlet(i).jacobianNotCalculated(); + } + // inlets + for (size_t i = 0; i < r->nInlets(); i++) { + r->inlet(i).jacobianNotCalculated(); + } + } + } // Create jacobian triplet vector vector jstarts; // Get jacobians and give elements to preconditioners @@ -600,6 +630,31 @@ void ReactorNet::buildJacobian(vector>& jacVector) jacVector[j] = newTrip; } } + + // loop through all connections and then set them found so calculations are not + // repeated + for (auto r : m_reactors) { + // walls + if (!m_jac_skip_walls) { + for (size_t i = 0; i < r->nWalls(); i++) { + r->wall(i).buildNetworkJacobian(jacVector); + r->wall(i).jacobianCalculated(); + } + } + // flow devices + if (m_jac_skip_flow_devices) { + // outlets + for (size_t i = 0; i < r->nOutlets(); i++) { + r->outlet(i).buildNetworkJacobian(jacVector); + r->outlet(i).jacobianCalculated(); + } + // inlets + for (size_t i = 0; i < r->nInlets(); i++) { + r->inlet(i).buildNetworkJacobian(jacVector); + r->inlet(i).jacobianCalculated(); + } + } + } } Eigen::SparseMatrix ReactorNet::finiteDifferenceJacobian() diff --git a/src/zeroD/Wall.cpp b/src/zeroD/Wall.cpp index b9287994847..e4ccaf30ccd 100644 --- a/src/zeroD/Wall.cpp +++ b/src/zeroD/Wall.cpp @@ -6,6 +6,8 @@ #include "cantera/base/stringUtils.h" #include "cantera/numerics/Func1.h" #include "cantera/zeroD/Wall.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/zeroD/ReactorNet.h" namespace Cantera { @@ -95,4 +97,66 @@ double Wall::heatRate() return q1; } + +void Wall::buildReactorJacobian(ReactorBase* r, vector>& jacVector) +{ + // get derivative of heat transfer for both reactors + vector> network; + size_t nsp = r->contents().nSpecies(); + size_t sidx = r->speciesOffset(); + size_t eidx = r->energyIndex(); + // define a scalar for direction based on left and right + double scalar = (r == m_left) ? 1.0 : -1.0; + // elements within the current reactor + // find dQdni for the current reactor w.r.t current reactor + for (size_t i = sidx; i < nsp + sidx; i++) { + double dQdni = m_rrth * m_area * scalar * r->moleDerivative(i); + dQdni += m_emiss * m_area * scalar * r->moleRadiationDerivative(i); + jacVector.emplace_back(eidx, i, dQdni); + } +} + +void Wall::buildNetworkJacobian(vector>& jacVector) +{ + // No interdependent terms for reservoirs + if (m_right->type() == "Reservoir" || m_left->type() == "Reservoir") { + return; + } + // return if the jacobian has been calculated + if (m_jac_calculated) { + return; + } + // get derivatives for inter-dependent reactor terms + //variables for the right side + vector> network; + size_t r_nsp = m_right->contents().nSpecies(); + size_t r_sidx = m_right->speciesOffset(); + size_t r_net = m_right->network().globalStartIndex((Reactor* ) m_right); + size_t r_eidx = m_right->energyIndex(); + + // variables for the left side + size_t l_nsp = m_left->contents().nSpecies(); + size_t l_sidx = m_left->speciesOffset(); + size_t l_net = m_left->network().globalStartIndex((Reactor* ) m_left); + size_t l_eidx = m_left->energyIndex(); + + if (((Reactor* ) m_right)->energyEnabled()) { + // find dQdni for the right reactor w.r.t left reactor + for (size_t i = l_sidx; i < l_sidx + l_nsp; i++) { + double dQdni = m_rrth * m_area * m_left->moleDerivative(i); + dQdni += m_emiss * m_area * m_left->moleRadiationDerivative(i); + jacVector.emplace_back(r_eidx + r_net, i + l_net, dQdni); + } + } + + if (((Reactor* ) m_left)->energyEnabled()) { + // find dQdni for the left reactor w.r.t right reactor + for (size_t i = r_sidx; i < r_sidx + r_nsp; i++) { + double dQdni = - m_rrth * m_area * m_right->moleDerivative(i); + dQdni -= m_emiss * m_area * m_right->moleRadiationDerivative(i); + jacVector.emplace_back(l_eidx + l_net, i + r_net, dQdni); + } + } +} + } diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index 5c119f1508f..b8ea200e89b 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -1195,21 +1195,67 @@ def test_get_solver_type(self): self.net2.initialize() self.assertEqual(self.net2.linear_solver_type, "GMRES") + def test_heat_transfer_network(self): + # create first reactor + gas1 = ct.Solution("h2o2.yaml", "ohmech") + gas1.TPX = 600, ct.one_atm, "O2:1.0" + r1 = self.reactorClass(gas1) + r1.chemistry_enabled = False + + # create second reactor + gas2 = ct.Solution("h2o2.yaml", "ohmech") + gas2.TPX = 300, ct.one_atm, "O2:1.0" + r2 = self.reactorClass(gas2) + r2.chemistry_enabled = False + + # create wall + U = 2.0 + A = 3.0 + w = ct.Wall(r1, r2, U=U, A=A) + net = ct.ReactorNet([r1, r2]) + jac = net.jacobian + + # check for values + for i in range(gas1.n_species): + assert np.isclose(jac[0, i + 1 + r1.n_vars], - U * A * gas2.T) + assert np.isclose(jac[r1.n_vars, i + 1], U * A * gas1.T) + if (i + 1 != 2): + assert np.isclose(jac[0, i + 1], U * A * gas1.T, rtol=1e-2) + assert np.isclose(jac[r1.n_vars, i + 1 + r1.n_vars], - U * A * gas2.T, rtol=1e-2) + class TestIdealGasMoleReactor(TestMoleReactor): reactorClass = ct.IdealGasMoleReactor test_preconditioner_unsupported = None - @pytest.mark.diagnose def test_heat_transfer_network(self): - # Result should be the same if (m * cp) / (U * A) is held constant - self.make_reactors(T1=300, T2=1000) - self.add_wall(U=200, A=1.0) - self.net.preconditioner = ct.AdaptivePreconditioner() - print(self.net.finite_difference_jacobian) - # self.net.advance(1.0) - # T1a = self.r1.T - # T2a = self.r2.T + # create first reactor + gas1 = ct.Solution("h2o2.yaml", "ohmech") + gas1.TPX = 600, ct.one_atm, "O2:1.0" + r1 = self.reactorClass(gas1) + r1.chemistry_enabled = False + + # create second reactor + gas2 = ct.Solution("h2o2.yaml", "ohmech") + gas2.TPX = 300, ct.one_atm, "O2:1.0" + r2 = self.reactorClass(gas2) + r2.chemistry_enabled = False + + # create wall + U = 2.0 + A = 3.0 + w = ct.Wall(r1, r2, U=U, A=A) + net = ct.ReactorNet([r1, r2]) + jac = net.jacobian + + # check for values + for i in range(gas1.n_species): + assert np.isclose(jac[0, i + 2 + r1.n_vars], - U * A * gas2.T) + assert np.isclose(jac[r1.n_vars, i + 2], U * A * gas1.T) + if (i + 2 != 3): + assert np.isclose(jac[0, i + 2], U * A * gas1.T, rtol=1e-2) + assert np.isclose(jac[r1.n_vars, i + 2 + r1.n_vars], - U * A * gas2.T, rtol=1e-2) + def test_adaptive_precon_integration(self): # Network one with non-mole reactor