diff --git a/data/transport/collisions.xml b/data/transport/collisions.xml index d86bf74e..44fd7296 100644 --- a/data/transport/collisions.xml +++ b/data/transport/collisions.xml @@ -2923,4 +2923,4 @@ 31.00 23.60 14.70 9.09 7.27 6.38 5.88 5.55 5.10 4.81 4.59 - \ No newline at end of file + diff --git a/src/general/Mixture.h b/src/general/Mixture.h index 76143b5e..9b9020f3 100644 --- a/src/general/Mixture.h +++ b/src/general/Mixture.h @@ -80,6 +80,22 @@ class Mixture state()->energyTransferSource(p_source); } + /** + * Provides density-based energy transfer source jacobian terms based on + * the current state of the mixture. + */ + void energyTransferJacobiansRho(double* const p_jacobianRho) { + state()->energyTransferJacobiansRho(p_jacobianRho); + } + + /** + * Provides temperature-based energy transfer source jacobian terms based on + * the current state of the mixture. + */ + void energyTransferJacobiansTTv(double* const p_jacobianTTv) { + state()->energyTransferJacobiansTTv(p_jacobianTTv); + } + /** * Add a named element composition to the mixture which may be retrieved * with getComposition(). diff --git a/src/kinetics/JacobianManager.cpp b/src/kinetics/JacobianManager.cpp index 4f718a0e..bb716414 100644 --- a/src/kinetics/JacobianManager.cpp +++ b/src/kinetics/JacobianManager.cpp @@ -79,6 +79,29 @@ void ReactionStoich::contributeToJacobian( //============================================================================== +template +void ReactionStoich::contributeToNetProgressJacobian( + const double kf, const double kb, const double* const conc, + double* const work, double* const p_dropRho, const size_t ns, const size_t r) const +{ + // Need to make sure that product species are zeroed out in the work + // array (don't need to zero out whole array) + for (int i = 0; i < ns; ++i) + work[i] = 0.0; + + // Compute the derivative of reaction rate with respect to the reactants + // and products (all other terms are zero) + m_reacs.diffRR(kf, conc, work, Equals()); + m_prods.diffRR(kb, conc, work, MinusEquals()); + + // Only loop over the necessary species + for (int i = 0; i < ns; ++i) + p_dropRho[i+r*ns] = work[i]; + +} + +//============================================================================== + template void ThirdbodyReactionStoich::contributeToJacobian( const double kf, const double kb, const double* const conc, @@ -104,6 +127,34 @@ void ThirdbodyReactionStoich::contributeToJacobian( //============================================================================== +template +void ThirdbodyReactionStoich::contributeToNetProgressJacobian( + const double kf, const double kb, const double* const conc, + double* const work, double* const p_dropRho, const size_t ns, const size_t r) const +{ + const double rrf = m_reacs.rr(kf, conc); + const double rrb = m_prods.rr(kb, conc); + const double rr = rrf - rrb; + double tb = 0.0; + + for (int i = 0; i < ns; ++i) + work[i] = 0.0; + + for (int i = 0; i < ns; ++i) { + work[i] = mp_alpha[i] * rr; + tb += mp_alpha[i] * conc[i]; + } + + m_reacs.diffRR(kf, conc, work, PlusEqualsTimes(tb)); + m_prods.diffRR(kb, conc, work, MinusEqualsTimes(tb)); + + for (int i = 0; i < ns; ++i) + p_dropRho[i+r*ns] = work[i]; + +} + +//============================================================================== + template void JacobianManager::addReactionStoich( JacStoichBase* p_reacs, JacStoichBase* p_prods, const StoichType type, @@ -293,6 +344,179 @@ void JacobianManager::addReaction(const Reaction& reaction) //============================================================================== +void JacobianManager::computeTReactionTDerivatives( + const std::vector& reactions, const double* const p_t, + std::vector >& dTrxndT) const + { + + for(size_t i = 0; i < reactions.size(); ++i) { + const size_t type = reactions[i].type(); + switch (type) { + case ASSOCIATIVE_IONIZATION: + dTrxndT.push_back({1.0, 0.0}); + break; + case DISSOCIATIVE_RECOMBINATION: + dTrxndT.push_back({0.0, 1.0}); + break; + case ASSOCIATIVE_DETACHMENT: + dTrxndT.push_back({1.0, 0.0}); + break; + case DISSOCIATIVE_ATTACHMENT: + dTrxndT.push_back({0.0, 1.0}); + break; + case DISSOCIATION_M: + dTrxndT.push_back({0.5*sqrt(p_t[0]*p_t[1])/p_t[0], 1.0}); + break; + case DISSOCIATION_E: + dTrxndT.push_back({0.0, 0.0}); + break; + case RECOMBINATION_M: + dTrxndT.push_back({1.0, 0.5*sqrt(p_t[0]*p_t[1])/p_t[0]}); + break; + case RECOMBINATION_E: + dTrxndT.push_back({0.0, 0.0}); + break; + case IONIZATION_M: + dTrxndT.push_back({1.0, 1.0}); + break; + case IONIZATION_E: + dTrxndT.push_back({0.0, 1.0}); + break; + case ION_RECOMBINATION_M: + dTrxndT.push_back({1.0, 1.0}); + break; + case ION_RECOMBINATION_E: + dTrxndT.push_back({1.0, 0.0}); + break; + case ELECTRONIC_ATTACHMENT_E: + dTrxndT.push_back({0.0, 0.0}); + break; + case ELECTRONIC_DETACHMENT_E: + dTrxndT.push_back({0.0, 0.0}); + break; + case ELECTRONIC_ATTACHMENT_M: + dTrxndT.push_back({0.0, 1.0}); + break; + case ELECTRONIC_DETACHMENT_M: + dTrxndT.push_back({1.0, 0.0}); + break; + case EXCHANGE: + dTrxndT.push_back({1.0, 1.0}); + break; + case EXCITATION_E: + dTrxndT.push_back({0.0, 0.0}); + break; + case EXCITATION_M: + dTrxndT.push_back({1.0, 1.0}); + break; + default: + dTrxndT.push_back({0.0, 0.0}); + break; + } + } + } + +//============================================================================== + +void JacobianManager::computeTReactionTvDerivatives( + const std::vector& reactions, const double* const p_t, + std::vector >& dTrxndTv) const + { + + for(size_t i = 0; i < reactions.size(); ++i) { + const size_t type = reactions[i].type(); + switch (type) { + case ASSOCIATIVE_IONIZATION: + dTrxndTv.push_back({0.0, 1.0}); + break; + case DISSOCIATIVE_RECOMBINATION: + dTrxndTv.push_back({1.0, 0.0}); + break; + case ASSOCIATIVE_DETACHMENT: + dTrxndTv.push_back({0.0, 1.0}); + break; + case DISSOCIATIVE_ATTACHMENT: + dTrxndTv.push_back({1.0, 0.0}); + break; + case DISSOCIATION_M: + dTrxndTv.push_back({0.5*sqrt(p_t[0]*p_t[1])/p_t[1], 0.0}); + break; + case DISSOCIATION_E: + dTrxndTv.push_back({1.0, 1.0}); + break; + case RECOMBINATION_M: + dTrxndTv.push_back({0.0, 0.5*sqrt(p_t[0]*p_t[1])/p_t[1]}); + break; + case RECOMBINATION_E: + dTrxndTv.push_back({1.0, 1.0}); + break; + case IONIZATION_M: + dTrxndTv.push_back({0.0, 0.0}); + break; + case IONIZATION_E: + dTrxndTv.push_back({1.0, 0.0}); + break; + case ION_RECOMBINATION_M: + dTrxndTv.push_back({0.0, 0.0}); + break; + case ION_RECOMBINATION_E: + dTrxndTv.push_back({0.0, 1.0}); + break; + case ELECTRONIC_ATTACHMENT_E: + dTrxndTv.push_back({1.0, 1.0}); + break; + case ELECTRONIC_DETACHMENT_E: + dTrxndTv.push_back({1.0, 1.0}); + break; + case ELECTRONIC_ATTACHMENT_M: + dTrxndTv.push_back({1.0, 0.0}); + break; + case ELECTRONIC_DETACHMENT_M: + dTrxndTv.push_back({0.0, 1.0}); + break; + case EXCHANGE: + dTrxndTv.push_back({0.0, 0.0}); + break; + case EXCITATION_E: + dTrxndTv.push_back({1.0, 1.0}); + break; + case EXCITATION_M: + dTrxndTv.push_back({0.0, 0.0}); + break; + default: + dTrxndTv.push_back({0.0, 0.0}); + break; + } + } + } + +//============================================================================== + +void JacobianManager::computeNetRateProgressJacobian( + const double* const kf, const double* const kb, const double* const conc, + double* const p_dropRho) const +{ + const size_t ns = m_thermo.nSpecies(); + const size_t nr = m_reactions.size(); + + // Make sure we are staring with a clean slate + std::fill(p_dropRho, p_dropRho + ns*nr, 0.0); + + // Loop over each reaction and compute the dRR/dconc_k + for (int i = 0; i < nr; ++i) + m_reactions[i]->contributeToNetProgressJacobian( + kf[i], kb[i], conc, mp_work1, p_dropRho, ns, i); + + + // Finally, multiply by the species molecular weight ratios + for (int i = 0, index = 0; i < nr; ++i) + for (int j = 0; j < ns; ++j, ++index) + p_dropRho[index] /= m_thermo.speciesMw(j); + +} + +//============================================================================== + void JacobianManager::computeJacobian( const double* const kf, const double* const kb, const double* const conc, double* const sjac) const @@ -304,17 +528,63 @@ void JacobianManager::computeJacobian( std::fill(sjac, sjac+ns*ns, 0.0); // Loop over each reaction and compute the dRR/dconc_k - for (int i = 0; i < nr; ++i) + for (int i = 0; i < nr; ++i) m_reactions[i]->contributeToJacobian( kf[i], kb[i], conc, mp_work, sjac, ns); + // Finally, multiply by the species molecular weight ratios - for (int i = 0, index = 0; i < ns; ++i) { + for (int i = 0, index = 0; i < ns; ++i) for (int j = 0; j < ns; ++j, ++index) sjac[index] *= m_thermo.speciesMw(i) / m_thermo.speciesMw(j); - } } +//============================================================================== + +void JacobianManager::computeJacobianT( + const double* const p_dropf, const double* const p_dropb, + const std::vector& reactions, double* const p_dropT) const + { + const size_t nt = m_thermo.nEnergyEqns(); + std::vector > dTrxndT; + + std::fill(mp_T, mp_T + nt, 0.); + m_thermo.getTemperatures(mp_T); + + if (nt == 2) { + computeTReactionTDerivatives(reactions, mp_T, dTrxndT); + } else { + for (int i = 0; i < reactions.size(); ++i) + dTrxndT.push_back({1.0, 1.0}); + } + + for (int i = 0; i < reactions.size(); ++i) + p_dropT[i] = p_dropf[i]*dTrxndT[i].first - p_dropb[i]*dTrxndT[i].second; + } + +//============================================================================== + +void JacobianManager::computeJacobianTv( + const double* const p_dropf, const double* const p_dropb, + const std::vector& reactions, double* const p_dropTv) const + { + const size_t nt = m_thermo.nEnergyEqns(); + std::vector > dTrxndTv; + + std::fill(mp_T, mp_T + nt, 0.); + m_thermo.getTemperatures(mp_T); + + if (nt == 2) { + computeTReactionTvDerivatives(reactions, mp_T, dTrxndTv); + } else { + for (int i = 0; i < reactions.size(); ++i) + dTrxndTv.push_back({0.0, 0.0}); + } + + for (int i = 0; i < reactions.size(); ++i) + p_dropTv[i] = p_dropf[i]*dTrxndTv[i].first - p_dropb[i]*dTrxndTv[i].second; + } + //============================================================================== } // namespace Kinetics diff --git a/src/kinetics/JacobianManager.h b/src/kinetics/JacobianManager.h index 1f8597a7..8365ddc2 100644 --- a/src/kinetics/JacobianManager.h +++ b/src/kinetics/JacobianManager.h @@ -329,6 +329,10 @@ class ReactionStoichBase virtual void contributeToJacobian( const double kf, const double kb, const double* const conc, double* const work, double* const sjac, const size_t ns) const = 0; + + virtual void contributeToNetProgressJacobian( + const double kf, const double kb, const double* const conc, + double* const work, double* const p_dropRho, const size_t ns, const size_t r) const = 0; }; /** @@ -379,6 +383,17 @@ class ReactionStoich : public ReactionStoichBase const double kf, const double kb, const double* const conc, double* const work, double* const sjac, const size_t ns) const; + /** + * Adds the net rate of progress Jacobian elements. Note + * that this will affect the rows of each species that particpates in this + * reaction. The compiler should be able to unroll each loop since the + * sizes are known at compile time given the Reactants and Products + * stoichiometry types. + */ + void contributeToNetProgressJacobian( + const double kf, const double kb, const double* const conc, + double* const work, double* const p_dropRho, const size_t ns, const size_t r) const; + protected: Reactants m_reacs; @@ -464,6 +479,16 @@ class ThirdbodyReactionStoich : public ReactionStoich const double kf, const double kb, const double* const conc, double* const work, double* const sjac, const size_t ns) const; + /** + * Computes the jacobians for net rates of progress. It is used to compute + * the jacobians for energy exchange due to ionization impact reactions. The + * size of the jacobians is nr*ns, where nr is number of reactions and + * ns is number of species. + */ + void contributeToNetProgressJacobian( + const double kf, const double kb, const double* const conc, + double* const work, double* const p_dropRho, const size_t ns, const size_t r) const; + friend void swap( ThirdbodyReactionStoich&, ThirdbodyReactionStoich&); @@ -491,6 +516,8 @@ class JacobianManager : m_thermo(thermo) { mp_work = new double [m_thermo.nSpecies()]; + mp_work1 = new double [m_thermo.nSpecies()]; + mp_T = new double [m_thermo.nEnergyEqns()]; } /** @@ -499,6 +526,8 @@ class JacobianManager ~JacobianManager() { delete [] mp_work; + delete [] mp_work1; + delete [] mp_T; std::vector::iterator iter = m_reactions.begin(); for ( ; iter != m_reactions.end(); ++iter) @@ -517,6 +546,45 @@ class JacobianManager void computeJacobian( const double* const kf, const double* const kb, const double* const conc, double* const sjac) const; + + /** + * Computes the net rate of progress Jacobian matrix. + */ + void computeNetRateProgressJacobian( + const double* const kf, const double* const kb, + const double* const conc, double* const npjac) const; + + /** + * Computes the translational-temperature-dependent + * chemistry Jacobian elements. + */ + void computeJacobianT( + const double* const p_ropf, const double* const p_ropb, + const std::vector& reactions, double* const p_dropT) const; + + /** + * Computes the vibrational-temperature-dependent + * chemistry Jacobian elements. + */ + void computeJacobianTv( + const double* const p_ropf, const double* const p_ropb, + const std::vector& reactions, double* const p_dropTv) const; + + /** + * Computes the reaction temperature derivatives + * with respect to translational temperature. + */ + void computeTReactionTDerivatives( + const std::vector& reactions, const double* const p_t, + std::vector >& dTrxndT) const; + + /** + * Computes the reaction temperature derivatives + * with respect to vibrational temperature. + */ + void computeTReactionTvDerivatives( + const std::vector& reactions, const double* const p_t, + std::vector >& dTrxndTv) const; private: @@ -544,6 +612,8 @@ class JacobianManager const Mutation::Thermodynamics::Thermodynamics& m_thermo; double* mp_work; + double* mp_work1; + double* mp_T; std::vector m_reactions; }; diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 9cdd55a1..b5d17b63 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -54,7 +54,12 @@ Kinetics::Kinetics( mp_ropf(NULL), mp_ropb(NULL), mp_rop(NULL), - mp_wdot(NULL) + mp_wdot(NULL), + mp_dropf(NULL), + mp_dropb(NULL), + mp_dropRho(NULL), + mp_dropT(NULL), + mp_dropTv(NULL) { if (mechanism == "none") return; @@ -104,6 +109,16 @@ Kinetics::~Kinetics() delete [] mp_rop; if (mp_wdot != NULL) delete [] mp_wdot; + if (mp_dropf != NULL) + delete [] mp_dropf; + if (mp_dropb != NULL) + delete [] mp_dropb; + if (mp_dropRho != NULL) + delete [] mp_dropRho; + if (mp_dropT != NULL) + delete [] mp_dropT; + if (mp_dropTv != NULL) + delete [] mp_dropTv; } //============================================================================== @@ -162,7 +177,11 @@ void Kinetics::closeReactions(const bool validate_mechanism) mp_ropb = new double [nReactions()]; mp_rop = new double [std::max(m_thermo.nSpecies(), (int) nReactions())]; mp_wdot = new double [m_thermo.nSpecies()]; - + mp_dropf = new double [nReactions()]; + mp_dropb = new double [nReactions()]; + mp_dropRho = new double [nReactions()*m_thermo.nSpecies()]; + mp_dropT = new double [nReactions()]; + mp_dropTv = new double [nReactions()]; } //============================================================================== @@ -279,6 +298,41 @@ void Kinetics::backwardRatesOfProgress( m_thirdbodies.multiplyThirdbodies(p_conc, p_ropb); } +//============================================================================== + +void Kinetics::forwardRateOfProgressDerivatives(double* const p_dropf) +{ + if(nReactions() == 0) + return; + + mp_rates->update(m_thermo); + ArrayXd dkfdT = + Map(mp_rates->dkfdT(), nReactions()); + + forwardRatesOfProgress(mp_ropf); + + for (int i = 0; i < nReactions(); ++i) + p_dropf[i] = mp_ropf[i]*dkfdT[i]; +} + +//============================================================================== + +void Kinetics::backwardRateOfProgressDerivatives(double* const p_dropb) +{ + if (nReactions() == 0) + return; + + mp_rates->update(m_thermo); + ArrayXd dkbdT = + Map(mp_rates->dkbdT(), nReactions()); + + backwardRatesOfProgress(mp_ropb); + + for(int i=0; i < nReactions(); ++i) + p_dropb[i] = mp_ropb[i]*dkbdT[i]; +} + + /* //============================================================================== @@ -367,6 +421,65 @@ void Kinetics::netProductionRates(double* const p_wdot) //============================================================================== +void Kinetics::netRateProgressRhoJacobian(double* const p_dropRho) +{ + // Special case of no reactions + if (nReactions() == 0) { + std::fill(p_dropRho, p_dropRho + m_thermo.nSpecies()*nReactions(), 0); + return; + } + + forwardRateCoefficients(mp_ropf); + backwardRateCoefficients(mp_ropb); + + // Compute species concentrations (mol/m^3) + Map(mp_rop, m_thermo.nSpecies()) = + (m_thermo.numberDensity() / NA) * + Map(m_thermo.X(), m_thermo.nSpecies()); + + // Compute the Jacobian matrix + m_jacobian.computeNetRateProgressJacobian(mp_ropf, mp_ropb, mp_rop, p_dropRho); +} + +//============================================================================== + +void Kinetics::netRateProgressTJacobian(double* const p_dropT) +{ + std::fill(p_dropT, p_dropT + nReactions(), 0); + + // Special case of no reactions + if (nReactions() == 0) return; + + + // Compute forward and backward rates of progress + forwardRateOfProgressDerivatives(mp_dropf); + backwardRateOfProgressDerivatives(mp_dropb); + + // Compute the Jacobian matrix + m_jacobian.computeJacobianT(mp_dropf, mp_dropb, m_reactions, p_dropT); + +} + +//============================================================================== + +void Kinetics::netRateProgressTvJacobian(double* const p_dropTv) +{ + std::fill(p_dropTv, p_dropTv + nReactions(), 0); + + // Special case of no reactions + if (nReactions() == 0) return; + + // Compute forward and backward rates of progress + forwardRateOfProgressDerivatives(mp_dropf); + backwardRateOfProgressDerivatives(mp_dropb); + + // Compute the Jacobian matrix + m_jacobian.computeJacobianTv(mp_dropf, mp_dropb, m_reactions, p_dropTv); + +} + +//============================================================================== + void Kinetics::jacobianRho(double* const p_jac) { // Special case of no reactions @@ -385,6 +498,58 @@ void Kinetics::jacobianRho(double* const p_jac) // Compute the Jacobian matrix m_jacobian.computeJacobian(mp_ropf, mp_ropb, mp_rop, p_jac); + +} + +//============================================================================== + +void Kinetics::jacobianT(double* const p_jacT) +{ + std::fill(p_jacT, p_jacT + m_thermo.nSpecies(), 0); + + // Special case of no reactions + if (nReactions() == 0) return; + + + // Compute forward and backward rates of progress + forwardRateOfProgressDerivatives(mp_dropf); + backwardRateOfProgressDerivatives(mp_dropb); + + // Compute the Jacobian matrix + m_jacobian.computeJacobianT(mp_dropf, mp_dropb, m_reactions, mp_dropT); + + m_reactants.decrSpecies(mp_dropT, p_jacT); + m_rev_prods.incrSpecies(mp_dropT, p_jacT); + m_irr_prods.incrSpecies(mp_dropT, p_jacT); + + for (int i = 0; i < m_thermo.nSpecies(); ++i) + p_jacT[i] *= m_thermo.speciesMw(i); + +} + +//============================================================================== + +void Kinetics::jacobianTv(double* const p_jacTv) +{ + std::fill(p_jacTv, p_jacTv + m_thermo.nSpecies(), 0); + + // Special case of no reactions + if (nReactions() == 0) return; + + // Compute forward and backward rates of progress + forwardRateOfProgressDerivatives(mp_dropf); + backwardRateOfProgressDerivatives(mp_dropb); + + // Compute the Jacobian matrix + m_jacobian.computeJacobianTv(mp_dropf, mp_dropb, m_reactions, mp_dropTv); + + m_reactants.decrSpecies(mp_dropTv, p_jacTv); + m_rev_prods.incrSpecies(mp_dropTv, p_jacTv); + m_irr_prods.incrSpecies(mp_dropTv, p_jacTv); + + for (int i = 0; i < m_thermo.nSpecies(); ++i) + p_jacTv[i] *= m_thermo.speciesMw(i); + } //============================================================================== diff --git a/src/kinetics/Kinetics.h b/src/kinetics/Kinetics.h index eb100139..7c07c0f0 100644 --- a/src/kinetics/Kinetics.h +++ b/src/kinetics/Kinetics.h @@ -149,6 +149,26 @@ class Kinetics void backwardRatesOfProgress( const double* const p_conc, double* const p_ropb); + /** + * Fills the vector dropf with the reaction temperature derivatives of + * forward rates of progress variables for each reaction + * \f$ \frac{dk_{f,j}}{dT_{reac}} \prod_i C_i^{\nu_{ij}^{'}} \Theta_{TB} \f$. + * + * @param p_dropf on return, the reaction temperature derivatives of the + * forward rates of progress in mol/m^3-s-K + */ + void forwardRateOfProgressDerivatives(double* const p_dropf); + + /** + * Fills the vector dropb with the reaction temperature derivatives of + * backward rates of progress variables for each reaction + * \f$ \frac{dk_{b,j}}{dT_{reac}} \prod_i C_i^{\nu_{ij}^{'}} \Theta_{TB} \f$. + * + * @param p_dropb on return, the reaction temperature derivatives of the + * backward rates of progress in mol/m^3-s-K + */ + void backwardRateOfProgressDerivatives(double* const p_dropb); + /** * Fills the vector rop with the net rates of progress for each reaction * \f$ \left[k_{f,j}\prod_i C_i^{\nu_{ij}^{'}}-k_{b,j}\prod_i @@ -181,6 +201,42 @@ class Kinetics */ void netProductionRates(double* const p_wdot); + /** + * Fills the matrix p_dropRho with the density-based net rate of progress jacobians. + * \f[ + * J_{i} = \frac{\partial \dot{\omega}_r}{\partial \rho_i} + * \f] + * The size of the Jacobian matrix is equal to the product of the number of species (ns) + * times the number of reactions (nr). + * + * @param p_dropRho - on return, the jacobian matrix \f$J_{i}\f$ + */ + void netRateProgressRhoJacobian(double* const p_dropRho); + + /** + * Fills the matrix p_dropT with the temperature-based net rate of + * progress jacobians. + * \f[ + * J_{i} = \frac{\partial \dot{\omega}_r}{\partial T_{tr}} + * \f] + * The size of the Jacobian matrix is equal to number of reactions. + * + * @param p_dropT - on return, the jacobian matrix \f$J_{i}\f$ + */ + void netRateProgressTJacobian(double* const p_dropT); + + /** + * Fills the matrix p_dropTv with the temperature-based net rate of + * progress jacobians. + * \f[ + * J_{i} = \frac{\partial \dot{\omega}_r}{\partial T_{ve}} + * \f] + * The size of the Jacobian matrix is equal to number of reactions. + * + * @param p_dropTv - on return, the jacobian matrix \f$J_{i}\f$ + */ + void netRateProgressTvJacobian(double* const p_dropTv); + /** * Fills the matrix p_jac with the species production rate jacobian matrix * \f[ @@ -194,6 +250,28 @@ class Kinetics */ void jacobianRho(double* const p_jac); + /** + * Fills the matrix p_jacT with the temperature jacobians for chemistry source term + * \f[ + * J_{i} = \frac{\partial \dot{\omega}_i}{\partial T_{tr}} + * \f] + * The Jacobian matrix size should be equal to number of species. + * + * @param p_jacT - on return, the jacobian matrix \f$J_{i}\f$ + */ + void jacobianT(double* const p_jacT); + + /** + * Fills the matrix p_jacTv with the temperature jacobians for chemistry source term + * \f[ + * J_{i} = \frac{\partial \dot{\omega}_i}{\partial T_{ve}} + * \f] + * The Jacobian matrix size should be equal to number of species. + * + * @param p_jacTv - on return, the jacobian matrix \f$J_{i}\f$ + */ + void jacobianTv(double* const p_jacTv); + /** * Returns the change in some species quantity across each reaction. */ @@ -237,6 +315,11 @@ class Kinetics double* mp_ropb; double* mp_rop; double* mp_wdot; + double* mp_dropf; + double* mp_dropb; + double* mp_dropRho; + double* mp_dropT; + double* mp_dropTv; }; diff --git a/src/kinetics/RateLawGroup.h b/src/kinetics/RateLawGroup.h index 78b7b679..33d51c3c 100644 --- a/src/kinetics/RateLawGroup.h +++ b/src/kinetics/RateLawGroup.h @@ -93,6 +93,12 @@ class RateLawGroup virtual void lnk( const Thermodynamics::StateModel* const p_state, double* const p_lnk) = 0; + /** + * Evaluates all of the temperature derivative of rates in the group and stores in the given vector. + */ + virtual void invkdkdT( + const Thermodynamics::StateModel* const p_state, double* const p_dkdT) = 0; + /** * Computes \Delta G / RT for this rate law group and subtracts these values * for each of the reactions in this group. @@ -101,14 +107,23 @@ class RateLawGroup { // Compute G_i/RT - ln(Patm/RT) const double val = std::log(ONEATM / (RU * m_t)); - for (int i = 0; i < ns; ++i) + for (size_t i = 0; i < ns; ++i) p_g[i] -= val; - + // Now subtract \Delta[G_i/RT - ln(Patm/RT)]_j m_reacs.decrReactions(p_g, p_r); m_prods.incrReactions(p_g, p_r); } + /** + * Computes 1 / Keq \frac{\partial Keq}{\partial T} for this rate law group + * and subtracts these values for each of the reactions in this group. + */ + void derivativeKeq(double* const p_dKeqdT, double* const p_dkdT) const + { + m_reacs.decrReactions(p_dKeqdT, p_dkdT); + m_prods.incrReactions(p_dKeqdT, p_dkdT); + } protected: @@ -172,6 +187,29 @@ class RateLawGroup1T : public RateLawGroup m_last_t = m_t; } + /** + * Evaluates all of the temperature derivative rates in the group and stores in the given vector. + */ + virtual void invkdkdT( + const Thermodynamics::StateModel* const p_state, double* const p_dkdT) + { + // Determine the reaction temperature for this group + m_t = TSelectorType().getT(p_state); + + // Update only if the temperature has changed + //if (std::abs(m_t - m_last_t) > 1.0e-10) { + const double invT = 1.0 / m_t; + + for (size_t i = 0; i < m_rates.size(); ++i) { + const std::pair& rate = m_rates[i]; + p_dkdT[rate.first] = rate.second.derivativebykf(invT); + } + //} + + // Save this temperature + m_last_t = m_t; + } + private: /// vector of rates to evaluate @@ -259,6 +297,21 @@ class RateLawGroupCollection iter->second->lnk(p_state, p_lnk); } + /** + * Computes the derivative of rate coefficients + * \f$ \frac{1}{k_{f,j}} \frac{dk_{f,j}}{dT_{reac}} \f$ + * in this collection and stores them in the vector at the index corresponding + * to their respective reaction. + */ + void derivativeOfRateCoefficients( + const Thermodynamics::StateModel* const p_state, double* const p_dkdT) + { + // Compute the forward rate constants + GroupMap::iterator iter = m_group_map.begin(); + for ( ; iter != m_group_map.end(); ++iter) + iter->second->invkdkdT(p_state, p_dkdT); + } + /** * Subtracts ln(keq) from the provided rate coefficients. */ @@ -275,6 +328,24 @@ class RateLawGroupCollection } } + /** + * Computes the temperature derivative of keq as follows: + * \f$ \frac{1}{keq} \frac{\partial keq}{\partial T} \f$ + */ + void derivativeKeq( + const Thermodynamics::Thermodynamics& thermo, double* const p_dKeqdT, double* const p_dkdT) + { + const size_t ns = thermo.nSpecies(); + GroupMap::iterator iter = m_group_map.begin(); + for ( ; iter != m_group_map.end(); ++iter) { + const RateLawGroup* p_group = iter->second; + thermo.speciesSTdGOverRT(p_group->getT(), p_dKeqdT); + for(size_t i = 0; i < ns; ++i) + p_dKeqdT[i] += 1./p_group->getT(); + p_group->derivativeKeq(p_dKeqdT, p_dkdT); + } + } + private: /// Collection of RateLawGroup objects diff --git a/src/kinetics/RateLaws.h b/src/kinetics/RateLaws.h index 8f9ac6b8..85a45458 100644 --- a/src/kinetics/RateLaws.h +++ b/src/kinetics/RateLaws.h @@ -78,6 +78,10 @@ class Arrhenius : public RateLaw return (k*invT*(m_n + m_temp*invT)); } + inline double derivativebykf(const double invT) const { + return (invT*(m_n + m_temp*invT)); + } + double A() const { return std::exp(m_lnA); } @@ -89,7 +93,7 @@ class Arrhenius : public RateLaw double T() const { return m_temp; } - + private: static std::vector sm_aunits; diff --git a/src/kinetics/RateManager.cpp b/src/kinetics/RateManager.cpp index 40d4460b..c778aa3c 100644 --- a/src/kinetics/RateManager.cpp +++ b/src/kinetics/RateManager.cpp @@ -45,6 +45,9 @@ public:\ inline double getT(const Thermodynamics::StateModel* const state) const {\ return ( __T__ );\ }\ + inline const char* getName(const Thermodynamics::StateModel* const state) const {\ + return #__NAME__;\ + }\ }; /// Temperature selector which returns the current translational temperature @@ -114,19 +117,22 @@ SELECT_RATE_LAWS(EXCITATION_E, ArrheniusTe, ArrheniusTe) RateManager::RateManager(size_t ns, const std::vector& reactions) : m_ns(ns), m_nr(reactions.size()), mp_lnkf(NULL), mp_lnkb(NULL), - mp_gibbs(NULL) + mp_gibbs(NULL), mp_dkfdT(NULL), mp_dkbdT(NULL), mp_dKeqdT(NULL) { // Add all of the reactions' rate coefficients to the manager const size_t nr = reactions.size(); for (size_t i = 0; i < m_nr; ++i) addReaction(i, reactions[i]); - + // Allocate storage in one block for both rate coefficient arrays and // species gibbs free energies - const size_t block_size = 2*m_nr + ns; + const size_t block_size = 4*m_nr + 2*ns; mp_lnkf = new double [block_size]; mp_lnkb = mp_lnkf + m_nr; mp_gibbs = mp_lnkb + m_nr; + mp_dkfdT = mp_gibbs + ns; + mp_dkbdT = mp_dkfdT + m_nr; + mp_dKeqdT = mp_dkbdT + m_nr; // Initialize the arrays to zero std::fill(mp_lnkf, mp_lnkf+block_size, 0.0); @@ -147,7 +153,7 @@ void RateManager::addReaction(const size_t rxn, const Reaction& reaction) { // Get the rate law which is being used in this reaction const RateLaw* p_rate = reaction.rateLaw(); - + // Arrhenius reactions if (typeid(*p_rate) == typeid(Arrhenius)) { selectRate(rxn, reaction); @@ -210,7 +216,7 @@ void RateManager::addRate(const size_t rxn, const Reaction& reaction) // note: mp_lnkff+(rxn+m_nr) = mp_lnkfb+rxn m_rate_groups.addRateCoefficient( rxn+m_nr, reaction.rateLaw()); - + m_rate_groups.addReaction(rxn, reaction); } else { @@ -225,16 +231,25 @@ void RateManager::update(const Thermodynamics::Thermodynamics& thermo) // Evaluate all of the different rate coefficients m_rate_groups.logOfRateCoefficients(thermo.state(), mp_lnkf); + // Evaluate all of the reaction temperature derivatives of + // rate coefficients + m_rate_groups.derivativeOfRateCoefficients(thermo.state(), mp_dkfdT); + // Copy rate coefficients which are the same as one of the previously // calculated ones std::vector::const_iterator iter = m_to_copy.begin(); for ( ; iter != m_to_copy.end(); ++iter) { const size_t index = *iter; mp_lnkb[index] = mp_lnkf[index]; + mp_dkbdT[index] = mp_dkfdT[index]; } // Subtract lnkeq(Tb) rate constants from the lnkf(Tb) to get lnkb(Tb) m_rate_groups.subtractLnKeq(thermo, mp_gibbs, mp_lnkb); + + // Evaluate the reaction temperature derivatives of equilibrium constants + m_rate_groups.derivativeKeq(thermo, mp_dKeqdT, mp_dkbdT); + } //============================================================================== diff --git a/src/kinetics/RateManager.h b/src/kinetics/RateManager.h index 5d93148c..365bd3d8 100644 --- a/src/kinetics/RateManager.h +++ b/src/kinetics/RateManager.h @@ -67,8 +67,8 @@ class RateManager const double* const lnkf() { return mp_lnkf; } /** - * Returns a pointer to the forward rate coefficients evaluated at the - * forward temperature. + * Returns a pointer to the backward rate coefficients evaluated at the + * backward temperature. */ const double* const lnkb() { return mp_lnkb; } @@ -79,6 +79,18 @@ class RateManager return m_irr; } + /** + * Returns a pointer to the derivatives of forward rate coefficients + * with respect to forward temperature. + */ + const double* const dkfdT() {return mp_dkfdT; } + + /** + * Returns a pointer to the derivatives of backward rate coefficients + * with respect to backward temperature. + */ + const double* const dkbdT() {return mp_dkbdT; } + private: /** @@ -122,6 +134,15 @@ class RateManager /// Storage for species Gibbs free energies double* mp_gibbs; + /// Storage for forward rate coefficient derivative at forward temperature + double* mp_dkfdT; + + /// Storage for forward rate coefficient derivative at backward temperature + double* mp_dkbdT; + + /// Storage for equilibrium constant temperature derivatives + double* mp_dKeqdT; + /// Stores the indices for which the forward and reverse temperature /// evaluations are equal std::vector m_to_copy; diff --git a/src/kinetics/StoichiometryManager.h b/src/kinetics/StoichiometryManager.h index bb9cf2c4..99366ece 100644 --- a/src/kinetics/StoichiometryManager.h +++ b/src/kinetics/StoichiometryManager.h @@ -110,7 +110,7 @@ class Stoich */ inline void decrReaction(const double* const p_s, double* const p_r) const { - for (int i = 0; i < N; ++i) + for (int i = 0; i < N; ++i) p_r[m_rxn] -= p_s[m_sps[i]]; } diff --git a/src/numerics/Functors.h b/src/numerics/Functors.h index 94354a6a..19b811c5 100644 --- a/src/numerics/Functors.h +++ b/src/numerics/Functors.h @@ -150,6 +150,19 @@ class PlusEqualsYDivAlpha const T m_a; }; +/// Equality functor x -= y / a +template +class MinusEqualsYDivAlpha +{ +public: + MinusEqualsYDivAlpha(const T& a) + : m_a(a) + {} + inline void operator()(T& x, const T& y) const { x -= y / m_a; } +private: + const T m_a; +}; + } // namespace Numerics } // namespace Mutation diff --git a/src/thermo/HarmonicOscillator.cpp b/src/thermo/HarmonicOscillator.cpp index 276394f5..de8294d7 100644 --- a/src/thermo/HarmonicOscillator.cpp +++ b/src/thermo/HarmonicOscillator.cpp @@ -33,13 +33,22 @@ namespace Mutation { double HarmonicOscillator::energy(double T) const { assert(T > 0.0); - double energy = 0.0; for (auto theta: m_characteristic_temps) energy += theta / (std::exp(theta/T) - 1.0); return energy; } +double HarmonicOscillator::energyTderivative(double T) const +{ + assert(T > 0.0); + double energyTderivative = 0.0; + for (auto theta: m_characteristic_temps) + energyTderivative += + std::pow(theta / T, 2) * std::exp(theta / T) / + std::pow(std::exp(theta / T) - 1.0, 2); + return energyTderivative; +} struct HarmonicOscillatorDB::Data { @@ -88,4 +97,4 @@ HarmonicOscillator HarmonicOscillatorDB::create( } } // namespace Thermodynamics -} // namespace Mutation \ No newline at end of file +} // namespace Mutation diff --git a/src/thermo/HarmonicOscillator.h b/src/thermo/HarmonicOscillator.h index 54870823..1fa9011f 100644 --- a/src/thermo/HarmonicOscillator.h +++ b/src/thermo/HarmonicOscillator.h @@ -45,6 +45,9 @@ class HarmonicOscillator /// Returns energy of vibrator in K at given temperature in K. double energy(double T) const; + /// Returns derivative of energy of vibrator with respect to temperature. + double energyTderivative(double T) const; + /// Convience function to access characteristic temperatures in K. const std::vector& characteristicTemperatures() const { return m_characteristic_temps; @@ -81,4 +84,4 @@ class HarmonicOscillatorDB } // namespace Mutation -#endif // THERMO_HARMONIC_OSCILLATOR_H \ No newline at end of file +#endif // THERMO_HARMONIC_OSCILLATOR_H diff --git a/src/thermo/MultiPhaseEquilSolver.h b/src/thermo/MultiPhaseEquilSolver.h index 4534d7b2..2651bcb0 100644 --- a/src/thermo/MultiPhaseEquilSolver.h +++ b/src/thermo/MultiPhaseEquilSolver.h @@ -32,7 +32,7 @@ #include #include - +using namespace Eigen; namespace Mutation { namespace Thermodynamics { diff --git a/src/thermo/Nasa7Polynomial.cpp b/src/thermo/Nasa7Polynomial.cpp index 6315e94a..5c6bc02d 100644 --- a/src/thermo/Nasa7Polynomial.cpp +++ b/src/thermo/Nasa7Polynomial.cpp @@ -104,6 +104,18 @@ void Nasa7Polynomial::gibbs(const double *const p_params, double &g) const mp_coefficients[index][6]; } +void Nasa7Polynomial::derivativeTgibbs(const double *const p_params, double &dg) const +{ + int index = tRange(p_params[0]); + + dg = mp_coefficients[index][0] * p_params[1] + + mp_coefficients[index][1] * p_params[2] + + mp_coefficients[index][2] * p_params[3] + + mp_coefficients[index][3] * p_params[4] + + mp_coefficients[index][4] * p_params[5] + + mp_coefficients[index][5] * p_params[6]; +} + void Nasa7Polynomial::computeParams(const double &T, double *const params, const ThermoFunction func) { @@ -146,6 +158,17 @@ void Nasa7Polynomial::computeParams(const double &T, double *const params, params[5] = -T4/20.0; params[6] = 1.0/T; } + + if (func == GIBBSPRIME) { + params[0] = T; + params[1] = -1.0/T; + params[2] = -0.5; + params[3] = -T/3.0; + params[4] = -T2*0.25; + params[5] = -T3*0.20; + params[6] = -1.0/T2; + } + } std::istream &operator>>(istream &in, Nasa7Polynomial &n7) diff --git a/src/thermo/Nasa7Polynomial.h b/src/thermo/Nasa7Polynomial.h index 4fd16f43..eadac6f3 100644 --- a/src/thermo/Nasa7Polynomial.h +++ b/src/thermo/Nasa7Polynomial.h @@ -147,6 +147,12 @@ class Nasa7Polynomial : public NasaPolynomial<2,7> */ void gibbs(const double *const p_params, double &g) const; + /** + * Computes temperature-derivative Gibbs free energy G/Ru/T. + * @see computeParams() + */ + void derivativeTgibbs(const double *const p_params, double &dg) const; + /** * Computes dimensionless specific heat, enthalpy, entropy, and gibbs free * energy - Cp/Ru, H/Ru/T, S/Ru, and G/Ru/T. @@ -164,7 +170,8 @@ class Nasa7Polynomial : public NasaPolynomial<2,7> CP, ENTHALPY, ENTROPY, - GIBBS + GIBBS, + GIBBSPRIME }; /** @@ -177,7 +184,8 @@ class Nasa7Polynomial : public NasaPolynomial<2,7> * CP - params[4] \n * ENTHALPY - params[6] \n * ENTROPY - params[5] \n - * GIBBS - params[7] + * GIBBS - params[7] \n + * GIBBSPRIME - params[7] * * Computing the parameters this way saves computation when they are reused * to compute thermodynamic properties for many species at once. diff --git a/src/thermo/Nasa9Polynomial.cpp b/src/thermo/Nasa9Polynomial.cpp index 5f52dd2e..c59d8d5e 100644 --- a/src/thermo/Nasa9Polynomial.cpp +++ b/src/thermo/Nasa9Polynomial.cpp @@ -112,10 +112,19 @@ void Nasa9Polynomial::gibbs(const double *const p_params, double &g) const int tr = tRange(-2.0 * p_params[3]); g = -mp_coefficients[tr][8]; - for (int i = 0; i < 8; ++i) + for (int i = 0; i < 8; ++i) g += mp_coefficients[tr][i] * p_params[i]; } +void Nasa9Polynomial::derivativeTgibbs(const double *const p_params, double &dg) const +{ + int tr = tRange(-3.0 * p_params[4]); + + dg = mp_coefficients[tr][0] * p_params[0]; + for (int i = 1; i < 8; ++i) + dg += mp_coefficients[tr][i] * p_params[i]; +} + void Nasa9Polynomial::computeParams(const double &T, double *const params, const ThermoFunction func) { @@ -172,6 +181,20 @@ void Nasa9Polynomial::computeParams(const double &T, double *const params, params[6] = -T4 / 20.0; params[7] = 1.0 / T; } + + // d(G(T)/RT)/dT = a1/(T^3) - a2*ln(T)/T - a3/T - a4/2 - a5*T/3 - a6*T^2/4 - + // a7*T^3/5 - b1/T^2 + if (func == GIBBSPRIME) { + params[0] = 1.0 / T3; + params[1] = -log(T) / T2; + params[2] = -1.0 / T; + params[3] = -0.5; + params[4] = -T / 3.0; + params[5] = -0.25 * T2; + params[6] = -T3 / 5.0; + params[7] = -1.0 / T2; + } + } int Nasa9Polynomial::tRange(double T) const diff --git a/src/thermo/Nasa9Polynomial.h b/src/thermo/Nasa9Polynomial.h index 82629ee8..51eb962e 100644 --- a/src/thermo/Nasa9Polynomial.h +++ b/src/thermo/Nasa9Polynomial.h @@ -95,6 +95,12 @@ class Nasa9Polynomial */ void gibbs(const double *const p_params, double &g) const; + /** + * Computes temperature derivative of Gibbs free energy G/Ru/T. + * @see computeParams() + */ + void derivativeTgibbs(const double *const p_params, double &g) const; + /** * Enumerates functions that require a parameters list. * @see computeParams() @@ -105,7 +111,8 @@ class Nasa9Polynomial ENTHALPY, ENTROPY, GIBBS, - THERMO + THERMO, + GIBBSPRIME }; /** @@ -119,7 +126,8 @@ class Nasa9Polynomial * ENTHALPY - params[8] \n * ENTROPY - params[7] \n * GIBBS - params[8] \n - * THERMO - params[22] + * THERMO - params[22] \n + * GIBBSPRIME - params[8] * * Computing the parameters this way saves computation when they are reused * to compute thermodynamic properties for many species at once. diff --git a/src/thermo/NasaDB.h b/src/thermo/NasaDB.h index fafaacf3..4d309d12 100644 --- a/src/thermo/NasaDB.h +++ b/src/thermo/NasaDB.h @@ -77,6 +77,11 @@ class NasaDB : public ThermoDB double* const g, double* const gt, double* const gr, double* const gv, double* const gel); + void derivativeTgibbs( + double Th, double Te, double Tr, double Tv, double Tel, double P, + double* const dg, double* const dgt, double* const dgr, double* const dgv, + double* const dgel); + protected: void loadAvailableSpecies(std::list& species_list); @@ -173,6 +178,24 @@ void NasaDB::gibbs( if (gel != NULL) std::fill(gel, gel+m_ns, 0.0); } +template +void NasaDB::derivativeTgibbs( + double Th, double Te, double Tr, double Tv, double Tel, double P, + double* const dg, double* const dgt, double* const dgr, double* const dgv, + double* const dgel) +{ + if (dg != NULL) { + PolynomialType::computeParams(Th, mp_params, PolynomialType::GIBBSPRIME); + for (size_t i = 0; i < m_ns; ++i) + m_polynomials[i].derivativeTgibbs(mp_params, dg[i]); + } + + if (dgt != NULL) std::fill(dgt, dgt+m_ns, 0.0); + if (dgr != NULL) std::fill(dgr, dgr+m_ns, 0.0); + if (dgv != NULL) std::fill(dgv, dgv+m_ns, 0.0); + if (dgel != NULL) std::fill(dgel, dgel+m_ns, 0.0); +} + template void NasaDB::loadAvailableSpecies( std::list& species_list) diff --git a/src/thermo/RrhoDB.cpp b/src/thermo/RrhoDB.cpp index 01d7cf4f..2e167886 100644 --- a/src/thermo/RrhoDB.cpp +++ b/src/thermo/RrhoDB.cpp @@ -447,6 +447,36 @@ class RrhoDB : public ThermoDB g[0] -= std::log(2.0); } + /** + * Computes the derivative of unitless Gibbs free energy of each species i, + * \f$\frac{\partial{G_i / R_u T_h}{\partial T_h}\f$ where \f$G_i\f$ is + * non-dimensionalized by the heavy particle translational temperature. + * + */ + void derivativeTgibbs( + double Th, double Te, double Tr, double Tv, double Tel, double P, + double* const dg, double* const dgt, double* const dgr, double* const dgv, + double* const dgel) + { + double Th2 = Th*Th; + + cpT(dg, EqDiv(Th)); + cpR(dg, PlusEqDiv(Th)); + cpV(Tv, dg, PlusEqDiv(Th)); + cpE(Tel, dg, PlusEqDiv(Th)); + + dsT(Th, Te, P, dg, MinusEq()); + dsR(Tr, dg, MinusEq()); + dsV(Tv, dg, MinusEq()); + dsE(Tel, dg, MinusEq()); + + hT(Th, Te, dg, MinusEqDiv(Th2)); + hR(Tr, dg, MinusEqDiv(Th2)); + hV(Tv, dg, MinusEqDiv(Th2)); + hE(Tel, dg, MinusEqDiv(Th2)); + + } + private: typedef Equals Eq; @@ -454,6 +484,7 @@ class RrhoDB : public ThermoDB typedef PlusEqualsYDivAlpha PlusEqDiv; typedef PlusEquals PlusEq; typedef MinusEquals MinusEq; + typedef MinusEqualsYDivAlpha MinusEqDiv; class ElecBFacsFunctor { @@ -813,11 +844,11 @@ class RrhoDB : public ThermoDB void sT(double Th, double Te, double P, double* const s, const OP& op) { double fac = 2.5 * (1.0 + std::log(Th)) - std::log(P); if (m_has_electron) - op(s[0], 2.5 * std::log(Te / Th) + fac + mp_lnqtmw[0]); + op(s[0], 2.5 * std::log(Te / Th) + fac + mp_lnqtmw[0]); for (int i = (m_has_electron ? 1 : 0); i < m_ns; ++i) op(s[i], fac + mp_lnqtmw[i]); } - + /** * Computes the unitless rotational entropy of each species. */ @@ -825,11 +856,11 @@ class RrhoDB : public ThermoDB void sR(double T, double* const s, const OP& op) { const double onelnT = 1.0 + std::log(T); LOOP_MOLECULES( - op(s[j], mp_rot_data[i].linearity * (onelnT - + op(s[j], mp_rot_data[i].linearity * (onelnT - mp_rot_data[i].ln_omega_t)); ) } - + /** * Computes the unitless vibrational entropy of each species. */ @@ -847,7 +878,7 @@ class RrhoDB : public ThermoDB op(s[j], (sum1 / T - sum2)); ) } - + /** * Computes the unitless electronic entropy of each species. */ @@ -866,6 +897,65 @@ class RrhoDB : public ThermoDB } } + /** + * Computes the temperature derivative of unitless translational entropy of each species. + */ + template + void dsT(double Th, double Te, double P, double* const ds, const OP& op) { + double fac = 2.5/Th; + if (m_has_electron) + op(ds[0], 0.0); + for (int i = (m_has_electron ? 1 : 0); i < m_ns; ++i) + op(ds[i], fac); + } + + /** + * Computes the temperature derivative of unitless rotational entropy of each species. + */ + template + void dsR(double T, double* const ds, const OP& op) { + const double onebyT = 1.0/T; + LOOP_MOLECULES( + op(ds[j], mp_rot_data[i].linearity * onebyT); + ) + } + + /** + * Computes the temperature derivative of unitless vibrational entropy of each species. + */ + template + void dsV(double T, double* const ds, const OP& op) { + int ilevel = 0; + double fac1, fac2, sum; + LOOP_MOLECULES( + sum = 0.0; + for (int k = 0; k < mp_nvib[i]; ++k, ilevel++) { + fac1 = std::exp(mp_vib_temps[ilevel] / T); + fac2 = (fac1-1.)*(fac1-1.); + sum += mp_vib_temps[ilevel]*mp_vib_temps[ilevel]*fac1/fac2; + } + op(ds[j], (sum / (T*T*T))); + ) + } + + /** + * Computes the temperature derivative of unitless electronic entropy of each species. + */ + template + void dsE(double T, double* const p_ds, const OP& op) { + updateElecBoltzmannFactors(T); + op(p_ds[0], 0.0); + + double* facs = mp_el_bfacs; + for (int i = 0; i < m_elec_data.nheavy; ++i, facs += 3) { + if (facs[0] > 0) + op(p_ds[i+m_elec_data.offset], + (facs[2]*facs[0]-facs[1]*facs[1])/(T*T*T*facs[0]*facs[0])); + else + op(p_ds[i+m_elec_data.offset], 0.0); + } + } + private: int m_ns; diff --git a/src/thermo/StateModel.h b/src/thermo/StateModel.h index 413ef0ae..d213cfb2 100644 --- a/src/thermo/StateModel.h +++ b/src/thermo/StateModel.h @@ -84,6 +84,8 @@ class StateModel mp_X = new double [m_thermo.nSpecies()]; for (int i = 0; i < thermo.nSpecies(); ++i) mp_X[i] = 0.0; + mp_omegaRho = new double [m_thermo.nSpecies()]; + mp_omegaTTv = new double [nenergy]; } /** @@ -93,6 +95,8 @@ class StateModel { // Delete data pointers delete [] mp_X; + delete [] mp_omegaRho; + delete [] mp_omegaTTv; // Delete transfer models for (int i = 0; i < m_transfer_models.size(); ++i) @@ -210,8 +214,8 @@ class StateModel /** * Returns a vector of length n_species times n_energies with each corresponding - * energy per unit mass. The first n_species values correspond to the total energy - * vector. The remaining n_species vectors correspond to the additional energy modes. + * energy per unit mass. The first n_species values correspond to the total energy + * vector. The remaining n_species vectors correspond to the additional energy modes. */ virtual void getEnergiesMass(double* const p_e) { throw NotImplementedError("StateModel::getEnergiesMass()"); @@ -219,8 +223,8 @@ class StateModel /** * Returns a vector of length n_species times n_energies with each corresponding - * enthalpy per unit mass. The first n_species values correspond to the total enthalpy - * vector. The remaining n_species vectors correspond to the additional energy modes. + * enthalpy per unit mass. The first n_species values correspond to the total enthalpy + * vector. The remaining n_species vectors correspond to the additional energy modes. */ virtual void getEnthalpiesMass(double* const p_h) { throw NotImplementedError("StateModel::getEnthalpiesMass()"); @@ -228,8 +232,8 @@ class StateModel /** * Returns a vector of length n_species times n_energies with each corresponding - * cp per unit mass. Each n_species vector corresponds to a temperature in the state - * model. The first one is associated with the heavy particle translational temperature. + * cp per unit mass. Each n_species vector corresponds to a temperature in the state + * model. The first one is associated with the heavy particle translational temperature. */ virtual void getCpsMass(double* const p_Cp){ @@ -238,16 +242,16 @@ class StateModel /** * Returns a vector of length n_species times n_energies with each corresponding - * cp per unit mass. Each n_species vector corresponds to a temperature in the state - * model. The first one is associated with the heavy particle translational temperature. + * cp per unit mass. Each n_species vector corresponds to a temperature in the state + * model. The first one is associated with the heavy particle translational temperature. */ virtual void getCvsMass(double* const p_Cp){ throw NotImplementedError("StateModel::getCvsMass()"); } - /** - * Assigns a unique temperature to each energy mode. + /** + * Assigns a unique temperature to each energy mode. */ virtual void getTagModes(int* const p_tag) { throw NotImplementedError("StateModel::getTagModes()"); @@ -280,7 +284,42 @@ class StateModel p_omega[m_transfer_models[i].first] += m_transfer_models[i].second->source(); } - + + /** + * This function provides the density-based jacobians for energy + * transfer source terms + */ + virtual void energyTransferJacobiansRho(double* const p_omegaJRho) + { + const int ns = m_thermo.nSpecies(); + std::fill(mp_omegaRho, mp_omegaRho + ns, 0.); + for (int i = 0; i < ns; ++i) + p_omegaJRho[i] = 0.0; + + for (int i = 0; i < m_transfer_models.size(); ++i) { + m_transfer_models[i].second->jacobianRho(mp_omegaRho); + for (int j = 0; j < ns; ++j) + p_omegaJRho[j] += mp_omegaRho[j]; + } + } + + /** + * This function provides the temperature-based jacobians for energy + * transfer source terms + */ + virtual void energyTransferJacobiansTTv(double* const p_omegaJTTv) + { + std::fill(mp_omegaTTv, mp_omegaTTv + m_nenergy, 0.); + for (int i = 0; i < m_nenergy; ++i) + p_omegaJTTv[i] = 0.0; + + for (int i = 0; i < m_transfer_models.size(); ++i) { + m_transfer_models[i].second->jacobianTTv(mp_omegaTTv); + for (int j = 0; j < m_nenergy; ++j) + p_omegaJTTv[j] += mp_omegaTTv[j]; + } + } + protected: /** * @todo add a removeTransferTerm for the case wehen the @@ -402,7 +441,8 @@ class StateModel double m_B; double* mp_X; - + double* mp_omegaRho; + double* mp_omegaTTv; std::vector< std::pair > m_transfer_models; diff --git a/src/thermo/ThermoDB.h b/src/thermo/ThermoDB.h index 37166863..926d1fdd 100644 --- a/src/thermo/ThermoDB.h +++ b/src/thermo/ThermoDB.h @@ -251,6 +251,28 @@ class ThermoDB double* const g, double* const gt, double* const gr, double* const gv, double* const gel) = 0; + /** + * Computes the temperature-derivative of Gibbs free energy of each species i, + * \f$\frac{\partial G_i / R_u T_h}{\partial T_h}\f$ where \f$G_i\f$ is + * non-dimensionalized by the heavy particle translational temperature. + * + * @param Th - heavy particle translational temperature + * @param Te - free electron temperature + * @param Tr - mixture rotational temperature + * @param Tv - mixture vibrational temperature + * @param Tel - mixture electronic temperature + * @param P - mixture static pressure + * @param dg - on return, the array of species energies derivatives + * @param dgt - if not NULL, the array of species translational energies derivatives + * @param dgr - if not NULL, the array of species rotational energies derivatives + * @param dgv - if not NULL, the array of species vibrational energies derivatives + * @param dgel - if not NULL, the array of species electronic energies derivatives + */ + virtual void derivativeTgibbs( + double Th, double Te, double Tr, double Tv, double Tel, double P, + double* const dg, double* const dgt, double* const dgr, double* const dgv, + double* const dgel) = 0; + protected: /** diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index ddad37ac..ee18df56 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -1045,6 +1045,13 @@ void Thermodynamics::speciesSTGOverRT(double T, double* const p_g) const { //============================================================================== +void Thermodynamics::speciesSTdGOverRT(double T, double* const p_dg) const { + mp_thermodb->derivativeTgibbs( + T, T, T, T, T, standardStateP(), p_dg, NULL, NULL, NULL, NULL); +} + +//============================================================================== + void Thermodynamics::elementMoles( const double *const species_N, double *const element_N) const { diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index 9ba1b888..8b6c862a 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -834,6 +834,13 @@ class Thermodynamics //: public StateModelUpdateHandler */ void speciesSTGOverRT(double T, double* const p_g) const; + /** + * Returns the temperature-derivative of species Gibbs free energies + * \f$ \frac{\partial G_i / R_u T}{\partial T} = + * \frac{\partial H_i / R_u T}{\partial T} - \frac{\partial S_i / R_u}{\partial T} \f$. + */ + void speciesSTdGOverRT(double T, double* const p_dg) const; + /** * Returns the number of moles of each element in a mixture with a given * set of species moles. diff --git a/src/transfer/OmegaCE.cpp b/src/transfer/OmegaCE.cpp index 7c516970..c40aca4a 100644 --- a/src/transfer/OmegaCE.cpp +++ b/src/transfer/OmegaCE.cpp @@ -40,26 +40,58 @@ namespace Mutation { class OmegaCE : public TransferModel { public: - OmegaCE(Mutation::Mixture& mix) - : TransferModel(mix) - { - mp_wrk1 = new double [mix.nSpecies()]; - } + OmegaCE(Mutation::Mixture& mix) + : TransferModel(mix) + { + mp_wrk1 = new double [mix.nSpecies()]; + mp_wrk2 = new double [mix.nSpecies()*mix.nSpecies()]; + mp_wrk3 = new double [mix.nSpecies()]; + mp_wrk4 = new double [mix.nSpecies()]; + } - ~OmegaCE() - { - delete [] mp_wrk1; - } + ~OmegaCE() + { + delete [] mp_wrk1; + delete [] mp_wrk2; + delete [] mp_wrk3; + delete [] mp_wrk4; + } - double source() - { - static const double cv = 1.5*RU/mixture().speciesMw(0); - mixture().netProductionRates(mp_wrk1); - return mp_wrk1[0]*cv*mixture().Te(); - } + double source() + { + static const double cv = 1.5*RU/mixture().speciesMw(0); + mixture().netProductionRates(mp_wrk1); + return mp_wrk1[0]*cv*mixture().Te(); + } + + void jacobianRho(double* const p_jacRho) + { + const size_t ns=mixture().nSpecies(); + std::fill(p_jacRho, p_jacRho + ns, 0.); + static const double cv = 1.5*RU/mixture().speciesMw(0); + mixture().jacobianRho(mp_wrk2); + for(int i = 0; i < ns; ++i) + p_jacRho[i] = mp_wrk2[i]*cv*mixture().Te(); + } + + void jacobianTTv(double* const p_jacTTv) + { + const size_t nt=mixture().nEnergyEqns(); + std::fill(p_jacTTv, p_jacTTv + nt, 0.); + static const double cv = 1.5*RU/mixture().speciesMw(0); + mixture().netProductionRates(mp_wrk1); + mixture().jacobianT(mp_wrk3); + mixture().jacobianTv(mp_wrk4); + p_jacTTv[0] = cv*mixture().Te()*mp_wrk3[0]; + p_jacTTv[1] = mp_wrk1[0]*cv; + p_jacTTv[1] += cv*mixture().Te()*mp_wrk4[0]; + } private: - double* mp_wrk1; + double* mp_wrk1; + double* mp_wrk2; + double* mp_wrk3; + double* mp_wrk4; }; // Register the transfer model diff --git a/src/transfer/OmegaCElec.cpp b/src/transfer/OmegaCElec.cpp index 84a11d76..1d20d9c1 100644 --- a/src/transfer/OmegaCElec.cpp +++ b/src/transfer/OmegaCElec.cpp @@ -43,35 +43,89 @@ namespace Mutation { class OmegaCElec : public TransferModel { public: - OmegaCElec(Mutation::Mixture& mix) - : TransferModel(mix) - { - mp_wrk1 = new double [mix.nSpecies()]; - mp_wrk2 = new double [mix.nSpecies()]; - }; - - ~OmegaCElec() - { - delete [] mp_wrk1; - delete [] mp_wrk2; - }; - - double source() - { - mixture().speciesHOverRT(NULL, NULL, NULL, NULL, mp_wrk1, NULL); - mixture().netProductionRates(mp_wrk2); - - double sum = 0.0; - - for(int i = 0 ; i < mixture().nSpecies(); ++i) - sum += mp_wrk1[i]*mp_wrk2[i]/mixture().speciesMw(i); - - return (sum*mixture().T()*RU); + OmegaCElec(Mutation::Mixture& mix) + : TransferModel(mix) + { + m_ns = mix.nSpecies(); + m_nt = mix.nEnergyEqns(); + mp_wrk1 = new double [mix.nSpecies()]; + mp_wrk2 = new double [mix.nSpecies()]; + mp_wrk3 = new double [mix.nSpecies()*mix.nSpecies()]; + mp_wrk4 = new double [mix.nSpecies()]; + mp_wrk5 = new double [mix.nSpecies()]; + mp_wrk6 = new double [mix.nSpecies()]; + }; + + ~OmegaCElec() + { + delete [] mp_wrk1; + delete [] mp_wrk2; + delete [] mp_wrk3; + delete [] mp_wrk4; + delete [] mp_wrk5; + delete [] mp_wrk6; + }; + + double source() + { + mixture().speciesHOverRT(NULL, NULL, NULL, NULL, mp_wrk1, NULL); + mixture().netProductionRates(mp_wrk2); + + double sum = 0.0; + + for(int i = 0 ; i < mixture().nSpecies(); ++i) + sum += mp_wrk1[i]*mp_wrk2[i]/mixture().speciesMw(i); + + return (sum*mixture().T()*RU); } + void jacobianRho(double* const p_jacRho) + { + std::fill(p_jacRho, p_jacRho + m_ns, 0.); + + mixture().speciesHOverRT(NULL, NULL, NULL, NULL, mp_wrk1, NULL); + mixture().jacobianRho(mp_wrk3); + + double aux = mixture().T()*RU; + + for(int i = 0; i < m_ns; ++i) + for(int j = 0; j < m_ns; ++j) + p_jacRho[i] += aux*mp_wrk1[i]*mp_wrk3[i+j*m_ns]/mixture().speciesMw(i); + + } + + void jacobianTTv(double* const p_jacTTv) + { + std::fill(p_jacTTv, p_jacTTv + m_nt, 0.); + + double T=mixture().T(); + double Tv=mixture().Tv(); + + mixture().speciesHOverRT(NULL, NULL, NULL, NULL, mp_wrk1, NULL); + mixture().netProductionRates(mp_wrk2); + mixture().speciesCpOverR(T, Tv, T, Tv, Tv, NULL, NULL, NULL, NULL, mp_wrk4); + mixture().jacobianT(mp_wrk5); + mixture().jacobianTv(mp_wrk6); + + double aux = T*RU; + + for(int i = 0; i < m_ns; ++i) { + p_jacTTv[0] += aux*mp_wrk1[i]*mp_wrk5[i]/mixture().speciesMw(i); + p_jacTTv[1] += aux*mp_wrk1[i]*mp_wrk6[i]/mixture().speciesMw(i); + p_jacTTv[1] += RU*mp_wrk4[i]*mp_wrk2[i]/mixture().speciesMw(i); + } + + } + private: - double* mp_wrk1; - double* mp_wrk2; + int m_ns; + int m_nt; + double* mp_wrk1; + double* mp_wrk2; + double* mp_wrk3; + double* mp_wrk4; + double* mp_wrk5; + double* mp_wrk6; }; // Register the transfer model diff --git a/src/transfer/OmegaCV.cpp b/src/transfer/OmegaCV.cpp index 4d3fdfa6..bd1089d8 100644 --- a/src/transfer/OmegaCV.cpp +++ b/src/transfer/OmegaCV.cpp @@ -56,44 +56,76 @@ class OmegaCV : public TransferModel */ public: - OmegaCV(Mutation::Mixture& mix) - : TransferModel(mix) - { - m_ns = mixture().nSpecies(); - mp_wrk1 = new double [m_ns]; - mp_wrk2 = new double [m_ns]; + OmegaCV(Mutation::Mixture& mix) + : TransferModel(mix) + { + m_ns = mixture().nSpecies(); + m_nt = mixture().nEnergyEqns(); + mp_wrk1 = new double [m_ns]; + mp_wrk2 = new double [m_ns]; + mp_wrk3 = new double [m_ns*m_ns]; + mp_wrk4 = new double [m_ns]; + mp_wrk5 = new double [m_ns]; + mp_wrk6 = new double [m_ns]; }; - ~OmegaCV() - { - delete [] mp_wrk1; - delete [] mp_wrk2; - }; + ~OmegaCV() + { + delete [] mp_wrk1; + delete [] mp_wrk2; + delete [] mp_wrk3; + delete [] mp_wrk4; + delete [] mp_wrk5; + delete [] mp_wrk6; + }; /** * Computes the source terms of the Vibration-Chemistry energy transfer in \f$ [J/(m^3\cdot s)] \f$ * * \f[ \Omega^{CV}_{mi} = \sum c_1 e^V_{mi} \dot{\omega}_i\f] * */ - double source() - { - static int i_transfer_model = 0; - switch (i_transfer_model){ - case 0: - return compute_source_Candler(); - break; - default: - std::cerr << "The selected Chemistry-Vibration-Chemistry model is not implemented yet"; - return 0.0; - } - } + double source() + { + static int i_transfer_model = 0; + switch (i_transfer_model){ + case 0: + return compute_source_Candler(); + break; + default: + std::cerr << "The selected Chemistry-Vibration-Chemistry model is not implemented yet"; + return 0.0; + } + } -private: - int m_ns; - double* mp_wrk1; - double* mp_wrk2; +/** + * Computes the density-based jacobian terms of the Vibration-Chemistry energy transfer.$ + * + * \f[ \frac{\partial \Omega^{CV}_{mi}}{\partial \rho_{i}} = + * \sum c_1 e^V_{mi} \frac{\partial \dot{\omega}_i}{\partial \rho_{i}} \f] + * + */ + void jacobianRho(double* const p_jacRho); + +/** + * Computes the temperature-based jacobian terms of the Vibration-Chemistry energy transfer.$ + * + * \f[ \frac{\partial \Omega^{CV}_{mi}}{\partial T} = + * \sum c_1 \frac{\partial e^V_{mi}}{\partial T} \dot{\omega}_i + * + \sum c_1 e^V_{mi} \frac{\partial \dot{\omega}_i}{\partial T} \f] + */ + void jacobianTTv(double* const p_jacTTv); - double const compute_source_Candler(); +private: + int m_ns; + int m_nt; + double* mp_wrk1; + double* mp_wrk2; + double* mp_wrk3; + double* mp_wrk4; + double* mp_wrk5; + double* mp_wrk6; + + double const compute_source_Candler(); }; /** @@ -106,25 +138,86 @@ class OmegaCV : public TransferModel * */ -double const OmegaCV::compute_source_Candler() -{ - + double const OmegaCV::compute_source_Candler() + { + + + // Getting Vibrational Energy + mixture().speciesHOverRT(NULL, NULL, NULL, mp_wrk1, NULL, NULL); + + // Getting Production Rate + mixture().netProductionRates(mp_wrk2); + + // Inner Product + double c1 = 1.0E0; + double sum = 0.E0; + + for(int i = 0 ; i < m_ns; ++i) + sum += mp_wrk1[i]*mp_wrk2[i]/mixture().speciesMw(i); + + return(c1*sum*mixture().T()*RU); + } - // Getting Vibrational Energy - mixture().speciesHOverRT(NULL, NULL, NULL, mp_wrk1, NULL, NULL); - - // Getting Production Rate - mixture().netProductionRates(mp_wrk2); +/** + * Assuming non-preferential model for now, density-based jacobians are computed. + * + */ - // Inner Product - double c1 = 1.0E0; - double sum = 0.E0; + void OmegaCV::jacobianRho(double* const p_jacRho) + { + std::fill(p_jacRho, p_jacRho + m_ns, 0.); + + // Getting Vibrational Energy + mixture().speciesHOverRT(NULL, NULL, NULL, mp_wrk1, NULL, NULL); + + // Getting density based jacobian for production rates + mixture().jacobianRho(mp_wrk3); + + //Inner product + double c1 = 1.0E0; + double aux = c1*RU*mixture().T(); + + for(int i = 0; i < m_ns; ++i) + for(int j = 0; j < m_ns; ++j) + p_jacRho[i] += aux*mp_wrk1[i]*mp_wrk3[j*m_ns+i]/mixture().speciesMw(i); + } - for(int i = 0 ; i < m_ns; ++i) - sum += mp_wrk1[i]*mp_wrk2[i]/mixture().speciesMw(i); +/** + * Assuming non-preferential model for now, temperature-based jacobians are computed. + * + */ - return(c1*sum*mixture().T()*RU); - } + void OmegaCV::jacobianTTv(double* const p_jacTTv) + { + std::fill(p_jacTTv, p_jacTTv + m_nt, 0.); + + double T=mixture().T(); + double Tv=mixture().Tv(); + + // Getting Vibrational Energy + mixture().speciesHOverRT(NULL, NULL, NULL, mp_wrk1, NULL, NULL); + + // Getting production rates + mixture().netProductionRates(mp_wrk2); + + // Getting vibrational energy derivative with respect to Tv + mixture().speciesCpOverR(T, Tv, T, Tv, Tv, NULL, NULL, NULL, mp_wrk4, NULL); + + // Getting temperature based jacobians for production rates + mixture().jacobianT(mp_wrk5); + mixture().jacobianTv(mp_wrk6); + + //Inner product + double c1 = 1.0E0; + double aux = c1*RU*T; + double aux1 = c1*RU; + + for(int i = 0; i < m_ns; ++i) { + p_jacTTv[0] += aux*mp_wrk1[i]*mp_wrk5[i]/mixture().speciesMw(i); + p_jacTTv[1] += aux*mp_wrk1[i]*mp_wrk6[i]/mixture().speciesMw(i); + p_jacTTv[1] += aux1*mp_wrk4[i]*mp_wrk2[i]/mixture().speciesMw(i); + } + } // Register the transfer model Mutation::Utilities::Config::ObjectProvider< diff --git a/src/transfer/OmegaET.cpp b/src/transfer/OmegaET.cpp index 71baf7ec..bcdebe68 100644 --- a/src/transfer/OmegaET.cpp +++ b/src/transfer/OmegaET.cpp @@ -57,38 +57,94 @@ namespace Mutation { class OmegaET : public TransferModel { public: - OmegaET(TransferModel::ARGS mix) : - TransferModel(mix), m_collisions(mix.collisionDB()), - m_has_electrons(mix.hasElectrons()) - { } - - double source() - { - if (!m_has_electrons) - return 0.0; - - const double * p_X = mixture().X(); + OmegaET(TransferModel::ARGS mix) : + TransferModel(mix), m_collisions(mix.collisionDB()), + m_has_electrons(mix.hasElectrons()) + { } + + double source() + { + if (!m_has_electrons) + return 0.0; + + const double * p_X = mixture().X(); double nd = mixture().numberDensity(); double T = mixture().T(); double Te = mixture().Te(); - + const double ns = mixture().nSpecies(); Map X(mixture().X(),ns); - + // CollisionDB data const ArrayXd& Q11ei = m_collisions.Q11ei(); const ArrayXd& mass = m_collisions.mass(); - + // Electron velocity double ve = sqrt(KB*8.*Te/(PI*mass(0))); double tau = 1.0/(mass(0)*8./3.*ve*nd*(X*Q11ei/mass).tail(ns-1).sum()); - + return 1.5*KB*nd*p_X[0]*(T-Te)/tau; - } - + } + + void jacobianRho(double* const p_jacRho) + { + std::fill(p_jacRho, p_jacRho + mixture().nSpecies(), 0.); + if (!m_has_electrons) return; + + const double * p_X = mixture().X(); + double nd = mixture().numberDensity(); + double T = mixture().T(); + double Te = mixture().Te(); + + const double ns = mixture().nSpecies(); + Map X(mixture().X(),ns); + + // CollisionDB data + const ArrayXd& Q11ei = m_collisions.Q11ei(); + const ArrayXd& mass = m_collisions.mass(); + + // Electron velocity + double ve = sqrt(KB*8.*Te/(PI*mass(0))); + double tau = 1.0/(mass(0)*8./3.*ve*nd*(X*Q11ei/mass).tail(ns-1).sum()); + + p_jacRho[0] = 1.5*KB*(NA/mixture().speciesMw(0))*(T-Te)/tau; + + double dtaudRho; + for(int i = 1; i < ns; ++i) { + dtaudRho = 1.0/(mass(0)*8./3.*ve*(NA/mixture().speciesMw(i))*Q11ei(i)/mass(i)); + p_jacRho[i] = 1.5*KB*nd*p_X[0]*(T-Te)/dtaudRho; + } + } + + void jacobianTTv(double* const p_jacTTv) + { + std::fill(p_jacTTv, p_jacTTv + mixture().nEnergyEqns(), 0.); + if (!m_has_electrons) return; + + const double* p_X = mixture().X(); + double nd = mixture().numberDensity(); + double T = mixture().T(); + double Te = mixture().Te(); + + const double ns = mixture().nSpecies(); + Map X(mixture().X(),ns); + + // CollisionDB data + const ArrayXd& Q11ei = m_collisions.Q11ei(); + const ArrayXd& mass = m_collisions.mass(); + + // Electron velocity + double ve = sqrt(KB*8.*Te/(PI*mass(0))); + double tau = 1.0/(mass(0)*8./3.*ve*nd*(X*Q11ei/mass).tail(ns-1).sum()); + + p_jacTTv[0] = 1.5*KB*nd*p_X[0]/tau; + p_jacTTv[1] = -1.5*KB*nd*p_X[0]/tau; + p_jacTTv[1] += (1.5*KB*nd*p_X[0]*(T-Te)/tau)*0.5/Te; + } + private: - Transport::CollisionDB& m_collisions; - bool m_has_electrons; + Transport::CollisionDB& m_collisions; + bool m_has_electrons; }; diff --git a/src/transfer/OmegaI.cpp b/src/transfer/OmegaI.cpp index 7b82b29a..ad95a74b 100644 --- a/src/transfer/OmegaI.cpp +++ b/src/transfer/OmegaI.cpp @@ -44,10 +44,16 @@ class OmegaI : public TransferModel { m_ns = mixture().nSpecies(); m_nr = mixture().nReactions(); + m_nt = mixture().nEnergyEqns(); mp_hf = new double [m_ns]; mp_h = new double [m_ns]; + mp_cp = new double [m_ns]; mp_rate = new double [m_nr]; mp_delta = new double [m_nr]; + mp_delta1 = new double [m_nr]; + mp_dropRho = new double [m_nr*m_ns]; + mp_dropT = new double [m_nr]; + mp_dropTv = new double [m_nr]; for(int i=0; i m_rId; double* mp_hf; double* mp_h; + double* mp_cp; double* mp_rate; double* mp_delta; + double* mp_delta1; + double* mp_dropRho; + double* mp_dropT; + double* mp_dropTv; }; // Register the transfer model diff --git a/src/transfer/OmegaVT.cpp b/src/transfer/OmegaVT.cpp index 11b7d04c..b429b240 100644 --- a/src/transfer/OmegaVT.cpp +++ b/src/transfer/OmegaVT.cpp @@ -62,6 +62,11 @@ class Vibrator return m_energy_model.energy(T); } + /// Temperature-derivative of energy of the vibrator at the given temperature. + double energyTderivative(double T) const { + return m_energy_model.energyTderivative(T); + } + /// Relaxation time of the vibrator. /// @todo Find better solution to limit dependency to just thermodynamic /// state information, not whole mixture. @@ -92,6 +97,8 @@ class OmegaVT : public TransferModel OmegaVT(Mixture&); virtual ~OmegaVT() = default; virtual double source() override; + virtual void jacobianRho(double* const p_jacRho) override; + virtual void jacobianTTv(double* const p_jacTTv) override; private: std::vector m_vibrators; }; @@ -138,9 +145,49 @@ double OmegaVT::source() return src * mixture().density() * RU; } +void OmegaVT::jacobianRho(double* const p_jacRho) +{ + const size_t ns=mixture().nSpecies(); + std::fill(p_jacRho, p_jacRho + ns, 0.); + auto Y = mixture().Y(); + auto T = mixture().T(); + auto Tv = mixture().Tv(); + for (auto& vibrator: m_vibrators) + { + const auto iv = vibrator.speciesIndex(); + const auto tau = vibrator.relaxationTime(mixture()); + const auto mw = vibrator.molecularWeight(); + p_jacRho[iv] = RU*(vibrator.energy(T) - vibrator.energy(Tv))/(mw*tau); + } +} + +void OmegaVT::jacobianTTv(double* const p_jacTTv) +{ + const size_t nt=mixture().nEnergyEqns(); + std::fill(p_jacTTv, p_jacTTv + nt, 0.); + + auto Y = mixture().Y(); + auto T = mixture().T(); + auto Tv = mixture().Tv(); + + double auxT = 0.0; + double auxTv = 0.0; + + for (auto& vibrator: m_vibrators) + { + const auto iv = vibrator.speciesIndex(); + const auto tau = vibrator.relaxationTime(mixture()); + const auto mw = vibrator.molecularWeight(); + auxT += Y[iv]*vibrator.energyTderivative(T)/(mw*tau); + auxTv += -Y[iv]*vibrator.energyTderivative(Tv)/(mw*tau); + } + + p_jacTTv[0] = auxT * mixture().density() * RU; + p_jacTTv[1] = auxTv * mixture().density() * RU; +} /** * Represents a coupling between vibrational and translational energy modes. diff --git a/src/transfer/TransferModel.h b/src/transfer/TransferModel.h index 67037a80..24d80f7f 100644 --- a/src/transfer/TransferModel.h +++ b/src/transfer/TransferModel.h @@ -61,6 +61,10 @@ class TransferModel /// Value of source term in J/m^3-s. virtual double source() = 0; + /// Value of jacobians for source term + virtual void jacobianRho(double* const p_jacRho) = 0; + virtual void jacobianTTv(double* const p_jacTTv) = 0; + protected: Mixture& mixture() { return m_mixture; } private: diff --git a/src/transport/CollisionDB.cpp b/src/transport/CollisionDB.cpp index aec89c86..7c082a50 100644 --- a/src/transport/CollisionDB.cpp +++ b/src/transport/CollisionDB.cpp @@ -243,6 +243,9 @@ const ArrayXd& CollisionDB::Dim(bool include_one_minus_x) const int k = ns - nh; const ArrayXd X = Map(m_thermo.X(), ns).max(1.0e-12); + // Single species - no diffusion by definition + if (ns == 1) { m_Dim.setZero(); return (m_Dim /= m_thermo.numberDensity()); } + // Electron if (k > 0) { nDei(); // update nDei @@ -251,6 +254,7 @@ const ArrayXd& CollisionDB::Dim(bool include_one_minus_x) } else m_Dim.setZero(); + // Heavies nDij(); // update nDij for (int i = 0, index = 1; i < nh; ++i, ++index) { diff --git a/tests/c++/test_reaction_mechanism.cpp b/tests/c++/test_reaction_mechanism.cpp index 38cfa152..14c2a7da 100644 --- a/tests/c++/test_reaction_mechanism.cpp +++ b/tests/c++/test_reaction_mechanism.cpp @@ -85,12 +85,20 @@ class TestMechanism Eigen::ArrayXd& forwardRatesOfProgress() { return m_rf; } /// Returns the reverse rates of progress. Eigen::ArrayXd& reverseRatesOfProgress() { return m_rb; } + /// Returns the reaction temperature derivatives of forward rates of progress. + Eigen::ArrayXd& forwardRatesOfProgressDerivatives() { return m_drf; } + /// Returns the reaction temperature derivatives of reverse rates of progress. + Eigen::ArrayXd& reverseRatesOfProgressDerivatives() { return m_drb; } /// Returns the reaction rates. Eigen::ArrayXd& netRatesOfProgress() { return m_rr; } /// Returns the production rates. Eigen::ArrayXd& productionRates() { return m_wdot; } /// Returns the production rate jacobians w.r.t. species densities. Eigen::MatrixXd& densityJacobian() { return m_jac_rho; } + /// Returns the production rate jacobians w.r.t translational temperature. + Eigen::ArrayXd& temperatureTrJacobian() { return m_jac_Ttr; } + /// Returns the production rate jacobians w.r.t vibrational temperature. + Eigen::ArrayXd& temperatureVeJacobian() { return m_jac_Tve; } private: @@ -102,11 +110,24 @@ class TestMechanism Eigen::ArrayXd m_kf; Eigen::ArrayXd m_kb; Eigen::ArrayXd m_keq; + Eigen::ArrayXd m_dkf; + Eigen::ArrayXd m_dkb; + Eigen::ArrayXd m_dkeq; + Eigen::ArrayXd m_dTfdTh; + Eigen::ArrayXd m_dTfdTv; + Eigen::ArrayXd m_dTbdTh; + Eigen::ArrayXd m_dTbdTv; Eigen::ArrayXd m_rf; Eigen::ArrayXd m_rb; + Eigen::ArrayXd m_drf; + Eigen::ArrayXd m_drb; + Eigen::ArrayXd m_drdTh; + Eigen::ArrayXd m_drdTv; Eigen::ArrayXd m_rr; Eigen::ArrayXd m_wdot; Eigen::MatrixXd m_jac_rho; + Eigen::ArrayXd m_jac_Ttr; + Eigen::ArrayXd m_jac_Tve; SharedPtr m_thermo; }; // class TestMechanism @@ -114,7 +135,9 @@ class TestMechanism TestMechanism::TestMechanism(bool reversible) : m_reversible(reversible), m_species("e- N O NO NO+ N2 O2"), m_A(5), m_Mw(7), m_kf(5), m_kb(5), m_keq(5), m_rf(5), m_rb(5), m_rr(5), - m_wdot(7), m_jac_rho(7, 7) + m_dkf(5), m_dkb(5), m_dkeq(5), m_dTfdTh(5), m_dTfdTv(5), m_dTbdTh(5), + m_dTbdTv(5), m_drf(5), m_drb(5), m_drdTh(5), m_drdTv(5), m_wdot(7), + m_jac_rho(7, 7), m_jac_Ttr(7), m_jac_Tve(7) { // Reaction formulas m_formulas.push_back("N2 + M = 2N + M"); @@ -179,6 +202,24 @@ void TestMechanism::setState(Array densities, double Th, double Tv) m_kf(2) = m_A(2) * Tv; m_kf(3) = m_A(3) * Th; m_kf(4) = m_A(4) * Th; + + m_dkf(0) = m_A(0); + m_dkf(1) = m_A(1); + m_dkf(2) = m_A(2); + m_dkf(3) = m_A(3); + m_dkf(4) = m_A(4); + + m_dTfdTh(0) = 0.5*std::sqrt(Th * Tv)/Th; + m_dTfdTh(1) = 0.5*std::sqrt(Th * Tv)/Th; + m_dTfdTh(2) = 0.; + m_dTfdTh(3) = 1.; + m_dTfdTh(4) = 1.; + + m_dTfdTv(0) = 0.5*std::sqrt(Th * Tv)/Tv; + m_dTfdTv(1) = 0.5*std::sqrt(Th * Tv)/Tv; + m_dTfdTv(2) = 1.; + m_dTfdTv(3) = 0.; + m_dTfdTv(4) = 0.; m_thermo->gibbs( Th, Th, Th, Th, Th, m_thermo->standardPressure(), @@ -192,16 +233,51 @@ void TestMechanism::setState(Array densities, double Th, double Tv) m_wdot.data(), NULL, NULL, NULL, NULL); m_keq(2) = ONEATM / (RU * Tv) * std::exp(-2.0*m_wdot(2) + m_wdot(6)); m_keq(4) = std::exp(-m_wdot(4) - m_wdot(0) + m_wdot(1) + m_wdot(2)); + + m_thermo->derivativeTgibbs( + Th, Th, Th, Th, Th, m_thermo->standardPressure(), + m_wdot.data(), NULL, NULL, NULL, NULL); + m_dkeq(0) = -1.0*m_keq(0)/Th - m_keq(0) * (2.0*m_wdot(1) - m_wdot(5)); + m_dkeq(1) = -1.0*m_keq(1)/Th - m_keq(1) * (2.0*m_wdot(2) - m_wdot(6)); + m_dkeq(3) = -m_keq(3) * (m_wdot(3) + m_wdot(2) - m_wdot(6) - m_wdot(1)); + + m_thermo->derivativeTgibbs( + Tv, Tv, Tv, Tv, Tv, m_thermo->standardPressure(), + m_wdot.data(), NULL, NULL, NULL, NULL); + m_dkeq(2) = -1.0*m_keq(2)/Tv - m_keq(2) * (2.0*m_wdot(2) - m_wdot(6)); + m_dkeq(4) = -m_keq(4) * (m_wdot(4) + m_wdot(0) - m_wdot(1) - m_wdot(2)); m_kb.setZero(); + m_dkb.setZero(); if (m_reversible) { m_kb(0) = m_A(0) * Th / m_keq(0); m_kb(1) = m_A(1) * Th / m_keq(1); m_kb(2) = m_A(2) * Tv / m_keq(2); m_kb(3) = m_A(3) * Th / m_keq(3); m_kb(4) = m_A(4) * Tv / m_keq(4); + m_dkb(0) = (m_A(0)*m_keq(0)-m_A(0)*Th*m_dkeq(0))/(m_keq(0)*m_keq(0)); + m_dkb(1) = (m_A(1)*m_keq(1)-m_A(1)*Th*m_dkeq(1))/(m_keq(1)*m_keq(1)); + m_dkb(2) = (m_A(2)*m_keq(2)-m_A(2)*Tv*m_dkeq(2))/(m_keq(2)*m_keq(2)); + m_dkb(3) = (m_A(3)*m_keq(3)-m_A(3)*Th*m_dkeq(3))/(m_keq(3)*m_keq(3)); + m_dkb(4) = (m_A(4)*m_keq(4)-m_A(4)*Tv*m_dkeq(4))/(m_keq(4)*m_keq(4)); } + m_dTbdTh.setZero(); + m_dTbdTv.setZero(); + if (m_reversible) { + m_dTbdTh(0) = 1.0; + m_dTbdTh(1) = 1.0; + m_dTbdTh(2) = 0.0; + m_dTbdTh(3) = 1.0; + m_dTbdTh(4) = 0.0; + + m_dTbdTv(0) = 0.0; + m_dTbdTv(1) = 0.0; + m_dTbdTv(2) = 1.0; + m_dTbdTv(3) = 0.0; + m_dTbdTv(4) = 1.0; + } + for (int i = 0; i < 7; ++i) m_wdot[i] = densities[i] / m_Mw[i]; double M = m_wdot.tail(6).sum(); @@ -212,6 +288,12 @@ void TestMechanism::setState(Array densities, double Th, double Tv) m_rf(3) = m_kf(3) * m_wdot(6) * m_wdot(1); m_rf(4) = m_kf(4) * m_wdot(1) * m_wdot(2); m_rr = m_rf; + + m_drf(0) = m_dkf(0) * m_wdot(5) * M; + m_drf(1) = m_dkf(1) * m_wdot(6) * m_wdot(1); + m_drf(2) = m_dkf(2) * m_wdot(6) * m_wdot(0); + m_drf(3) = m_dkf(3) * m_wdot(6) * m_wdot(1); + m_drf(4) = m_dkf(4) * m_wdot(1) * m_wdot(2); m_rb(0) = m_kb(0) * m_wdot(1) * m_wdot(1) * M; m_rb(1) = m_kb(1) * m_wdot(2) * m_wdot(2) * m_wdot(1); @@ -220,6 +302,15 @@ void TestMechanism::setState(Array densities, double Th, double Tv) m_rb(4) = m_kb(4) * m_wdot(4) * m_wdot(0); if (m_reversible) m_rr -= m_rb; + m_drb(0) = m_dkb(0) * m_wdot(1) * m_wdot(1) * M; + m_drb(1) = m_dkb(1) * m_wdot(2) * m_wdot(2) * m_wdot(1); + m_drb(2) = m_dkb(2) * m_wdot(2) * m_wdot(2) * m_wdot(0); + m_drb(3) = m_dkb(3) * m_wdot(3) * m_wdot(2); + m_drb(4) = m_dkb(4) * m_wdot(4) * m_wdot(0); + + m_drdTh = m_drf*m_dTfdTh - m_drb*m_dTbdTh; + m_drdTv = m_drf*m_dTfdTv - m_drb*m_dTbdTv; + // Compute jacobian terms Eigen::VectorXd drr(7), stoich(7); m_jac_rho.setConstant(0.0); @@ -280,6 +371,27 @@ void TestMechanism::setState(Array densities, double Th, double Tv) for (int i = 0; i < 7; ++i) for (int j = 0; j < 7; ++j) m_jac_rho(i,j) *= (m_Mw(i)/m_Mw(j)); + + m_jac_Ttr(0) = m_drdTh(4); + m_jac_Ttr(1) = 2*m_drdTh(0) - m_drdTh(3) - m_drdTh(4); + m_jac_Ttr(2) = 2*m_drdTh(1) + 2*m_drdTh(2) + m_drdTh(3) - m_drdTh(4); + m_jac_Ttr(3) = m_drdTh(3); + m_jac_Ttr(4) = m_drdTh(4); + m_jac_Ttr(5) = -m_drdTh(0); + m_jac_Ttr(6) = -m_drdTh(1) - m_drdTh(2) - m_drdTh(3); + + m_jac_Ttr *= m_Mw; + + m_jac_Tve(0) = m_drdTv(4); + m_jac_Tve(1) = 2*m_drdTv(0) - m_drdTv(3) - m_drdTv(4); + m_jac_Tve(2) = 2*m_drdTv(1) + 2*m_drdTv(2) + m_drdTv(3) - m_drdTv(4); + m_jac_Tve(3) = m_drdTv(3); + m_jac_Tve(4) = m_drdTv(4); + m_jac_Tve(5) = -m_drdTv(0); + m_jac_Tve(6) = -m_drdTv(1) - m_drdTv(2) - m_drdTv(3); + + m_jac_Tve *= m_Mw; + } /// Jumps to another permutation of x. @@ -328,6 +440,12 @@ void testMech(bool isRev) mix->backwardRatesOfProgress(test_nr.data()); REQUIRE(test_nr.isApprox(mech.reverseRatesOfProgress())); + + mix->forwardRateOfProgressDerivatives(test_nr.data()); + REQUIRE(test_nr.isApprox(mech.forwardRatesOfProgressDerivatives())); + + mix->backwardRateOfProgressDerivatives(test_nr.data()); + REQUIRE(test_nr.isApprox(mech.reverseRatesOfProgressDerivatives())); mix->netRatesOfProgress(test_nr.data()); REQUIRE(test_nr.isApprox(mech.netRatesOfProgress())); @@ -337,6 +455,12 @@ void testMech(bool isRev) mix->jacobianRho(test_jac.data()); REQUIRE(test_jac.isApprox(mech.densityJacobian())); + + mix->jacobianT(test_ns.data()); + REQUIRE(test_ns.isApprox(mech.temperatureTrJacobian())); + + mix->jacobianTv(test_ns.data()); + REQUIRE(test_ns.isApprox(mech.temperatureVeJacobian())); } } } while (stillPermuting(rho)); diff --git a/tests/c++/test_reaction_rates.cpp b/tests/c++/test_reaction_rates.cpp index 727df6e4..a25b2976 100644 --- a/tests/c++/test_reaction_rates.cpp +++ b/tests/c++/test_reaction_rates.cpp @@ -67,12 +67,12 @@ TEST_CASE("Arrhenius rate laws work correctly", "[kinetics]") { // Check that the rate is calculated correctly for (int order = 1; order < 4; ++order) { - checkArrheniusRate(5.0e15, 0.0, order, 75500.0); - checkArrheniusRate(5.0e15, 0.0, order, 75500.0, true); - checkArrheniusRate(5.0e15, 1.0, order, 75500.0); - checkArrheniusRate(5.0e15, 1.0, order, 75500.0, true); - checkArrheniusRate(5.0e15, -1.0, order, 75500.0); - checkArrheniusRate(5.0e15, -1.0, order, 75500.0, true); + checkArrheniusRate(5.0e15, 0.0, 75500.0, order); + checkArrheniusRate(5.0e15, 0.0, 75500.0, order, true); + checkArrheniusRate(5.0e15, 1.0, 75500.0, order); + checkArrheniusRate(5.0e15, 1.0, 75500.0, order, true); + checkArrheniusRate(5.0e15, -1.0, 75500.0, order); + checkArrheniusRate(5.0e15, -1.0, 75500.0, order, true); } // Check that loading fails properly