From 3193fa3906f2f90b5f6fb6c7fe78effabbfa4d35 Mon Sep 17 00:00:00 2001 From: Anthony Walker Date: Sun, 7 Apr 2024 10:27:46 -0400 Subject: [PATCH] Temp PR updates --- doc/sphinx/develop/index.md | 23 +++++------ include/cantera/zeroD/FlowDevice.h | 12 +++--- .../zeroD/IdealGasConstPressureMoleReactor.h | 4 +- include/cantera/zeroD/IdealGasMoleReactor.h | 4 +- include/cantera/zeroD/Reactor.h | 13 ------- include/cantera/zeroD/ReactorBase.h | 38 +++++++++++-------- include/cantera/zeroD/ReactorNet.h | 2 +- include/cantera/zeroD/Wall.h | 32 +++++++--------- .../IdealGasConstPressureMoleReactor.cpp | 18 ++------- src/zeroD/IdealGasMoleReactor.cpp | 22 ++--------- src/zeroD/Reactor.cpp | 2 +- src/zeroD/ReactorNet.cpp | 16 ++++---- src/zeroD/Wall.cpp | 21 +++++----- test/zeroD/test_zeroD.cpp | 8 ++-- 14 files changed, 88 insertions(+), 127 deletions(-) diff --git a/doc/sphinx/develop/index.md b/doc/sphinx/develop/index.md index 75daf963035..f30863de1ce 100644 --- a/doc/sphinx/develop/index.md +++ b/doc/sphinx/develop/index.md @@ -1,6 +1,7 @@ # Develop (sec-compiling)= + ## Compiling Cantera from Source If you're interested in contributing new features to Cantera, or you want to try the @@ -12,9 +13,9 @@ Cantera](compiling/configure-build) on your computer. The following additional references may also be useful: -- [](compiling/dependencies.md) -- [](compiling/config-options) -- [](compiling/special-cases) +- [](compiling/dependencies.md) +- [](compiling/config-options) +- [](compiling/special-cases) ```{toctree} :caption: Compiling Cantera from Source @@ -35,7 +36,7 @@ compiling/special-cases This section is a work in progress. ``` -- [](reactor-integration) +- [](reactor-integration) ```{toctree} :caption: How Cantera Works @@ -47,13 +48,13 @@ reactor-integration ## Adding New Features to Cantera -- [](CONTRIBUTING) -- [](style-guidelines) -- [](vscode-tips) -- [](writing-tests) -- [](running-tests) -- [](writing-examples) -- [](doc-formatting) +- [](CONTRIBUTING) +- [](style-guidelines) +- [](vscode-tips) +- [](writing-tests) +- [](running-tests) +- [](writing-examples) +- [](doc-formatting) ```{toctree} :caption: Adding New Features to Cantera diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index ffa61f66b71..36c1794c784 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -123,7 +123,7 @@ class FlowDevice //! @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. + //! @since New in %Cantera 3.1. //! virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { throw NotImplementedError(type() + "::buildReactorJacobian"); @@ -137,7 +137,7 @@ class FlowDevice //! @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. + //! @since New in %Cantera 3.1. //! virtual void buildNetworkJacobian(vector>& jacVector) { if (!m_jac_calculated) { @@ -149,16 +149,16 @@ class FlowDevice //! @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. + //! @since New in %Cantera 3.1. //! - void jacobianCalculated() { m_jac_calculated = true; }; + void calculatedJacobian() { 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. + //! @since New in %Cantera 3.1. //! - void jacobianNotCalculated() { m_jac_calculated = false; }; + void notCalculatedJacobian() { m_jac_calculated = false; }; protected: //! a variable to switch on and off so calculations are not doubled by the calling diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index d54060060c0..db186fc7bac 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -47,9 +47,7 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor bool preconditionerSupported() const override { return true; }; - double moleDerivative(size_t index) override; - - double moleRadiationDerivative(size_t index) override; + double temperature_ddni(size_t index) override; size_t speciesOffset() const override { return m_sidx; }; diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index fe245d84077..c5de1081320 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -43,9 +43,7 @@ class IdealGasMoleReactor : public MoleReactor bool preconditionerSupported() const override {return true;}; - double moleDerivative(size_t index) override; - - double moleRadiationDerivative(size_t index) override; + double temperature_ddni(size_t index) override; size_t speciesOffset() const override { return m_sidx; }; diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index dd433a42b30..37975d08377 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -84,12 +84,6 @@ class Reactor : public ReactorBase m_chem = cflag; } - //! Returns `true` if changes in the reactor composition due to chemical reactions - //! are enabled. - bool chemistryEnabled() const { - return m_chem; - } - void setEnergy(int eflag=1) override { if (eflag > 0) { m_energy = true; @@ -98,11 +92,6 @@ class Reactor : public ReactorBase } } - //! Returns `true` if solution of the energy equation is enabled. - bool energyEnabled() const { - return m_energy; - } - //! Number of equations (state variables) for this reactor size_t neq() { if (!m_nv) { @@ -335,8 +324,6 @@ class Reactor : public ReactorBase 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 diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 76c95c982ee..c0f6884c1c4 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -277,26 +277,15 @@ class ReactorBase //! Set the ReactorNet that this reactor belongs to. void setNetwork(ReactorNet* net); - //! Calculate the derivative of T with respect to the ith species in the heat - //! transfer equation based on the reactor specific equation of state. + //! Calculate the derivative of T with respect to the ith species in the energy + //! conservation 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. + //! @since New in %Cantera 3.1. //! - virtual double moleRadiationDerivative(size_t index) { - throw NotImplementedError("Reactor::moleRadiationDerivative"); + virtual double temperature_ddni(size_t index) { + throw NotImplementedError("Reactor::temperature_ddni"); } //! Return the index associated with energy of the system @@ -305,6 +294,17 @@ class ReactorBase //! Return the offset between species and state variables virtual size_t speciesOffset() const { return m_sidx; }; + //! Returns `true` if solution of the energy equation is enabled. + virtual bool energyEnabled() const { + return m_energy; + } + + //! Returns `true` if changes in the reactor composition due to chemical reactions + //! are enabled. + bool chemistryEnabled() const { + return m_chem; + } + protected: //! Specify the mixture contained in the reactor. Note that a pointer to //! this substance is stored, and as the integration proceeds, the state of @@ -348,6 +348,12 @@ class ReactorBase //! Composite thermo/kinetics/transport handler shared_ptr m_solution; + + //! A bool that enables the energy equation + bool m_energy = true; + + //! A bool that enables the chemical kinetics equations + bool m_chem = false; }; } diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 2b295a85a95..5db8665de9b 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -239,7 +239,7 @@ class ReactorNet : public FuncEval //! 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) { + size_t globalStartIndex(ReactorBase* curr_reactor) { for (size_t i = 0; i < m_reactors.size(); i++) { if (curr_reactor == m_reactors[i]) { return m_start[i]; diff --git a/include/cantera/zeroD/Wall.h b/include/cantera/zeroD/Wall.h index 548f75faaa8..d6ad965d3dd 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -95,12 +95,11 @@ class WallBase //! 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 + //! @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. + //! changed or removed without notice. + //! @since New in %Cantera 3.1. //! virtual void buildReactorJacobian(ReactorBase* r, vector>& jacVector) { throw NotImplementedError("WallBase::buildReactorJacobian"); @@ -109,32 +108,29 @@ class WallBase //! 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 + //! @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. + //! changed or removed without notice. + //! @since New in %Cantera 3.1. //! virtual void buildNetworkJacobian(vector>& jacVector) { throw NotImplementedError("WallBase::buildNetworkJacobian"); } - //! Specify the jacobian terms have been calculated and should not be recalculated. + //! 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. + //! changed or removed without notice. + //! @since New in %Cantera 3.1. //! - void jacobianCalculated() { m_jac_calculated = true; }; + void calculatedJacobian() { m_jac_calculated = true; }; - //! Specify that jacobian terms have not been calculated and should be recalculated. + //! 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. + //! changed or removed without notice. + //! @since New in %Cantera 3.1. //! - void jacobianNotCalculated() { m_jac_calculated = false; }; + void notCalculatedJacobian() { m_jac_calculated = false; }; protected: ReactorBase* m_left = nullptr; diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index f02c99f2496..b80cb01fe5a 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -273,23 +273,11 @@ string IdealGasConstPressureMoleReactor::componentName(size_t k) { "Index is out of bounds."); } -double IdealGasConstPressureMoleReactor::moleDerivative(size_t index) +double IdealGasConstPressureMoleReactor::temperature_ddni(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; + double n_total = m_mass / m_thermo->meanMolecularWeight(); + return pressure() * m_vol / GasConstant / n_total; } } diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index 98728778b93..452129d7558 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -244,32 +244,16 @@ 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) +double IdealGasMoleReactor::temperature_ddni(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; + double n_total = m_mass / m_thermo->meanMolecularWeight(); + return pressure() * m_vol / GasConstant / n_total; } } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index da2b7f6d70d..53db141c772 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -603,7 +603,7 @@ void Reactor::setAdvanceLimit(const string& nm, const double limit) } } -void Reactor:: buildWallJacobian(vector>& jacVector) +void Reactor::buildWallJacobian(vector>& jacVector) { if (!m_jac_skip_walls) { for (size_t i = 0; i < m_wall.size(); i++) { diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index dd4cb18de82..bb2ed8ae0f3 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -590,7 +590,7 @@ void ReactorNet::checkPreconditionerSupported() const { void ReactorNet::buildJacobian(vector>& jacVector) { // network must be initialized for the jacobian - if (! m_init) { + if (!m_init) { initialize(); } // loop through and set connectors not found @@ -598,18 +598,18 @@ void ReactorNet::buildJacobian(vector>& jacVector) // walls if (!m_jac_skip_walls) { for (size_t i = 0; i < r->nWalls(); i++) { - r->wall(i).jacobianNotCalculated(); + r->wall(i).notCalculatedJacobian(); } } // flow devices if (!m_jac_skip_flow_devices) { // outlets for (size_t i = 0; i < r->nOutlets(); i++) { - r->outlet(i).jacobianNotCalculated(); + r->outlet(i).notCalculatedJacobian(); } // inlets for (size_t i = 0; i < r->nInlets(); i++) { - r->inlet(i).jacobianNotCalculated(); + r->inlet(i).notCalculatedJacobian(); } } } @@ -638,7 +638,7 @@ void ReactorNet::buildJacobian(vector>& jacVector) if (!m_jac_skip_walls) { for (size_t i = 0; i < r->nWalls(); i++) { r->wall(i).buildNetworkJacobian(jacVector); - r->wall(i).jacobianCalculated(); + r->wall(i).calculatedJacobian(); } } // flow devices @@ -646,12 +646,12 @@ void ReactorNet::buildJacobian(vector>& jacVector) // outlets for (size_t i = 0; i < r->nOutlets(); i++) { r->outlet(i).buildNetworkJacobian(jacVector); - r->outlet(i).jacobianCalculated(); + r->outlet(i).calculatedJacobian(); } // inlets for (size_t i = 0; i < r->nInlets(); i++) { r->inlet(i).buildNetworkJacobian(jacVector); - r->inlet(i).jacobianCalculated(); + r->inlet(i).calculatedJacobian(); } } } @@ -664,7 +664,7 @@ Eigen::SparseMatrix ReactorNet::finiteDifferenceJacobian() initialize(); } - // clear former jacobian elements + // allocate jacobian triplet vector vector> jac_trips; // Get the current state diff --git a/src/zeroD/Wall.cpp b/src/zeroD/Wall.cpp index d24b2cb27ac..742dab58701 100644 --- a/src/zeroD/Wall.cpp +++ b/src/zeroD/Wall.cpp @@ -86,12 +86,13 @@ void Wall::buildReactorJacobian(ReactorBase* r, vector>& 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; + double direction = (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); + double dQdni = m_rrth * m_area * direction * r->temperature_ddni(i); + dQdni += m_emiss * m_area * direction * r->temperature_ddni(i) * 4 + * pow(r->temperature(), 3); jacVector.emplace_back(eidx, i, dQdni); } } @@ -111,20 +112,21 @@ void Wall::buildNetworkJacobian(vector>& jacVector) 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_net = m_right->network().globalStartIndex(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_net = m_left->network().globalStartIndex(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); + double dQdni = m_rrth * m_area * m_left->temperature_ddni(i); + dQdni += m_emiss * m_area * m_left->temperature_ddni(i) * 4 + * pow(m_left->temperature(), 3); jacVector.emplace_back(r_eidx + r_net, i + l_net, dQdni); } } @@ -132,8 +134,9 @@ void Wall::buildNetworkJacobian(vector>& jacVector) 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); + double dQdni = - m_rrth * m_area * m_right->temperature_ddni(i); + dQdni -= m_emiss * m_area * m_right->temperature_ddni(i) * 4 + * pow(m_right->temperature(), 3); jacVector.emplace_back(l_eidx + l_net, i + r_net, dQdni); } } diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index 44a1a97ef73..4aef1be230a 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -267,10 +267,10 @@ TEST(JacobianTests, test_wall_jacobian_build) } // check that switch works wallJac.clear(); - w.jacobianCalculated(); + w.calculatedJacobian(); w.buildNetworkJacobian(wallJac); EXPECT_EQ(wallJac.size(), 0); - w.jacobianNotCalculated(); + w.notCalculatedJacobian(); // build jac for network terms wallJac.clear(); w.buildNetworkJacobian(wallJac); @@ -319,10 +319,10 @@ TEST(JacobianTests, test_flow_jacobian_build) // expect errors from building jacobians EXPECT_THROW(mfc.buildReactorJacobian(&reactor, flowJac), NotImplementedError); // check the jacobian calculated flag and throw/catch errors accordingly - mfc.jacobianCalculated(); + mfc.calculatedJacobian(); mfc.buildNetworkJacobian(flowJac); EXPECT_EQ(flowJac.size(), 0); - mfc.jacobianNotCalculated(); + mfc.notCalculatedJacobian(); EXPECT_THROW(mfc.buildNetworkJacobian(flowJac), NotImplementedError); }