From 87e1c1bbe7c18ea6f3a11dce8158308e433a7489 Mon Sep 17 00:00:00 2001 From: Kyle Hanquist Date: Mon, 20 Jan 2025 14:01:00 -0700 Subject: [PATCH 01/23] New file from Catarina that helps with NAN in SU2 coupling, See 1/20/25 email --- src/thermo/ChemNonEqTTvStateModel.cpp | 46 ++++++++++----------------- 1 file changed, 17 insertions(+), 29 deletions(-) diff --git a/src/thermo/ChemNonEqTTvStateModel.cpp b/src/thermo/ChemNonEqTTvStateModel.cpp index 3a30acc3..c94e2b68 100644 --- a/src/thermo/ChemNonEqTTvStateModel.cpp +++ b/src/thermo/ChemNonEqTTvStateModel.cpp @@ -74,22 +74,19 @@ class ChemNonEqTTvStateModel : public StateModel * following: * 0: {species densities}, {total energy density, vib. energy density} * 1: {species densities} and {T, Tv} - * 2: {species mass fractions} and {P, T, Tv} array */ void setState( const double* const p_mass, const double* const p_energy, - const int vars = 0) + const int vars = 0, const double Tv_old = 0, const double T_old = 0) { const int ns = m_thermo.nSpecies(); // Compute the species concentrations which are used through out this - // method regardless of variable set. Note that in the case where - // the p_mass vector is mass fractions, this is just dividing by the - // molecular weights as the first step in converting to mole fractions. + // method regardless of variable set double conc = 0.0; for (int i = 0; i < ns; ++i){ // Check that species densities are at least positive - assert(p_mass[i] >= 0.0 && "Species is negative!"); + assert(p_mass[i] >= 0.0); mp_X[i] = std::max(p_mass[i] / m_thermo.speciesMw(i), 0.0); conc += mp_X[i]; } @@ -101,50 +98,38 @@ class ChemNonEqTTvStateModel : public StateModel switch (vars) { case 0: { - solveEnergies(p_mass, p_energy); + solveEnergies(p_mass, p_energy, Tv_old, T_old); break; } case 1: // Check that temperature is at least positive - assert(p_energy[0] > 0.0 && "T is negative!"); - assert(p_energy[1] > 0.0 && "Tv is negative!"); + assert(p_energy[0] > 0.0); + assert(p_energy[1] > 0.0); m_T = p_energy[0]; m_Tv = p_energy[1]; break; - case 2: - // Check temperature and pressure are positive - assert(p_energy[0] > 0.0 && "P is negative!"); - assert(p_energy[1] > 0.0 && "T is negative!"); - assert(p_energy[2] > 0.0 && "Tv is negative!"); - m_P = p_energy[0]; - m_T = p_energy[1]; - m_Tv = p_energy[2]; - break; default: throw InvalidInputError("variable set", vars) - << "This variable-set is not implemented in ChemNonEqTTv" + << "This variable-set is not implemented in ChemNonEqStateModel" << ". Possible variable-sets are:\n" << " 0: (species densities, static energy density)\n" - << " 1: (species densities, T and Tv)\n" - << " 2: (species mass fractions, P, T and Tv)"; + << " 1: (species densities, T and Tv)"; } // Set Tr, Tel, and Te m_Tr = m_T; m_Tel = m_Te = m_Tv; - // Compute the species mole fractions + // Compute the pressure and species mole fractions from T and rho_i for (int i = 0; i < ns; ++i) mp_X[i] /= conc; - - // Compute the pressure from T and rho_i - if (vars != 2) - m_P = conc * RU * (m_T + (m_Tv - m_T) * elec / conc); + m_P = conc * RU * (m_T + (m_Tv - m_T) * elec / conc); } - /** - * Initialize transfer models + * Add proper + * 0: {species densities}, {total energy density, vib. energy density} + * 1: {species densities} and {T, Tv} */ void initializeTransferModel(Mutation::Mixture& mix) { @@ -254,7 +239,7 @@ class ChemNonEqTTvStateModel : public StateModel * @param p_rhoe - total and internal mass energy */ - void solveEnergies(const double* const p_rhoi, const double* const p_rhoe) + void solveEnergies(const double* const p_rhoi, const double* const p_rhoe, const double Tv_old = 0, const double T_old = 0) { // This line is required to have perfect numerical jacobians... //m_T = m_Tv = 500.0; @@ -295,6 +280,9 @@ class ChemNonEqTTvStateModel : public StateModel f = e - emix; } + if (m_T != m_T) m_T = T_old; + if (m_Tv != m_Tv) m_Tv = Tv_old; + if (i == imax) std::cout << "Warning, didn't converge temperatures: f = " << f.norm() << std::endl; } From e436e0c4ec0ad29cd14da93bdc54b40ed90fb007 Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Mon, 20 Jan 2025 17:30:02 -0700 Subject: [PATCH 02/23] Successfully builds; added function arguments to setState() --- src/thermo/ChemNonEqStateModel.cpp | 2 +- src/thermo/EquilStateModel.cpp | 2 +- src/thermo/StateModel.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/thermo/ChemNonEqStateModel.cpp b/src/thermo/ChemNonEqStateModel.cpp index ad04d109..949bc5cf 100644 --- a/src/thermo/ChemNonEqStateModel.cpp +++ b/src/thermo/ChemNonEqStateModel.cpp @@ -62,7 +62,7 @@ class ChemNonEqStateModel : public StateModel */ void setState( const double* const p_mass, const double* const p_energy, - const int vars = 0) + const int vars = 0, const double Tv_old = 0, const double T_old = 0) { const int ns = m_thermo.nSpecies(); diff --git a/src/thermo/EquilStateModel.cpp b/src/thermo/EquilStateModel.cpp index c8155cb9..ae8fca1a 100644 --- a/src/thermo/EquilStateModel.cpp +++ b/src/thermo/EquilStateModel.cpp @@ -76,7 +76,7 @@ class EquilStateModel : public StateModel */ virtual void setState( const double* const p_mass, const double* const p_energy, - const int vars = 0) + const int vars = 0, const double Tv_old = 0, const double T_old = 0) { // Determine temperature, pressure, and mole fractions depending on // variable set diff --git a/src/thermo/StateModel.h b/src/thermo/StateModel.h index 413ef0ae..89bcd15a 100644 --- a/src/thermo/StateModel.h +++ b/src/thermo/StateModel.h @@ -120,7 +120,7 @@ class StateModel */ virtual void setState( const double* const p_mass, const double* const p_energy, - const int vars = 0) = 0; + const int vars = 0, const double Tv_old = 0, const double T_old = 0) = 0; /** * Sets the current magnitude of the magnetic field in teslas. From 382b64dd55577ab45d49070ece6409733427adf9 Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Tue, 21 Jan 2025 08:59:31 -0700 Subject: [PATCH 03/23] Updated function headers to include Tv_old and T_old to protect against NAN in SU2 --- src/thermo/EquilStateModel.cpp | 2 +- src/thermo/Thermodynamics.cpp | 2 +- src/thermo/Thermodynamics.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/thermo/EquilStateModel.cpp b/src/thermo/EquilStateModel.cpp index ae8fca1a..3c069009 100644 --- a/src/thermo/EquilStateModel.cpp +++ b/src/thermo/EquilStateModel.cpp @@ -272,7 +272,7 @@ class EquilTPStateModel : public EquilStateModel * Sets the state of the equilibrium mixture given temperature and pressure. */ virtual void setState( - const double* const p_T, const double* const p_P, const int vars = 1) + const double* const p_T, const double* const p_P, const int vars = 1, const double Tv_old = 0, const double T_old = 0) { //assert(vars == 1); EquilStateModel::setState(p_P, p_T, 1); diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index ddad37ac..5396e9e0 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -176,7 +176,7 @@ bool Thermodynamics::speciesThermoValidAtT(const size_t i, const double T) const //============================================================================== void Thermodynamics::setState( - const double* const p_v1, const double* const p_v2, const int vars) + const double* const p_v1, const double* const p_v2, const int vars, const double Tv_old = 0, const double T_old = 0) { mp_state->setState(p_v1, p_v2, vars); convert(X(), mp_y); diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index 9ba1b888..aad3c9fa 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -306,7 +306,7 @@ class Thermodynamics //: public StateModelUpdateHandler * used. */ void setState( - const double* const p_v1, const double* const p_v2, const int vars = 0); + const double* const p_v1, const double* const p_v2, const int vars = 0, const double Tv_old = 0, const double T_old = 0); /** * Computes the equilibrium composition of the mixture at the given fixed From bd2f4824386866a39e7b853618e227a843e799dc Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Tue, 21 Jan 2025 09:02:09 -0700 Subject: [PATCH 04/23] Fixed redefinition of default arguments --- src/thermo/Thermodynamics.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index 5396e9e0..32106175 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -176,7 +176,7 @@ bool Thermodynamics::speciesThermoValidAtT(const size_t i, const double T) const //============================================================================== void Thermodynamics::setState( - const double* const p_v1, const double* const p_v2, const int vars, const double Tv_old = 0, const double T_old = 0) + const double* const p_v1, const double* const p_v2, const int vars, const double Tv_old, const double T_old) { mp_state->setState(p_v1, p_v2, vars); convert(X(), mp_y); From 2b41465161e0570c8cba562ca36482a52911644a Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Tue, 21 Jan 2025 10:40:35 -0700 Subject: [PATCH 05/23] Builds successfully. Removes T_old, Tv_old --- src/thermo/Thermodynamics.cpp | 2 +- src/thermo/Thermodynamics.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index 32106175..ddad37ac 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -176,7 +176,7 @@ bool Thermodynamics::speciesThermoValidAtT(const size_t i, const double T) const //============================================================================== void Thermodynamics::setState( - const double* const p_v1, const double* const p_v2, const int vars, const double Tv_old, const double T_old) + const double* const p_v1, const double* const p_v2, const int vars) { mp_state->setState(p_v1, p_v2, vars); convert(X(), mp_y); diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index aad3c9fa..9ba1b888 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -306,7 +306,7 @@ class Thermodynamics //: public StateModelUpdateHandler * used. */ void setState( - const double* const p_v1, const double* const p_v2, const int vars = 0, const double Tv_old = 0, const double T_old = 0); + const double* const p_v1, const double* const p_v2, const int vars = 0); /** * Computes the equilibrium composition of the mixture at the given fixed From e8e719dfd7c9262cc04c9c26d8f22e9f4c983785 Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Mon, 3 Feb 2025 16:00:48 -0700 Subject: [PATCH 06/23] Revert "Successfully builds; added function arguments to setState()" This reverts commit e436e0c4ec0ad29cd14da93bdc54b40ed90fb007. --- src/thermo/ChemNonEqStateModel.cpp | 2 +- src/thermo/EquilStateModel.cpp | 2 +- src/thermo/StateModel.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/thermo/ChemNonEqStateModel.cpp b/src/thermo/ChemNonEqStateModel.cpp index 949bc5cf..ad04d109 100644 --- a/src/thermo/ChemNonEqStateModel.cpp +++ b/src/thermo/ChemNonEqStateModel.cpp @@ -62,7 +62,7 @@ class ChemNonEqStateModel : public StateModel */ void setState( const double* const p_mass, const double* const p_energy, - const int vars = 0, const double Tv_old = 0, const double T_old = 0) + const int vars = 0) { const int ns = m_thermo.nSpecies(); diff --git a/src/thermo/EquilStateModel.cpp b/src/thermo/EquilStateModel.cpp index 3c069009..cecda860 100644 --- a/src/thermo/EquilStateModel.cpp +++ b/src/thermo/EquilStateModel.cpp @@ -76,7 +76,7 @@ class EquilStateModel : public StateModel */ virtual void setState( const double* const p_mass, const double* const p_energy, - const int vars = 0, const double Tv_old = 0, const double T_old = 0) + const int vars = 0) { // Determine temperature, pressure, and mole fractions depending on // variable set diff --git a/src/thermo/StateModel.h b/src/thermo/StateModel.h index 89bcd15a..413ef0ae 100644 --- a/src/thermo/StateModel.h +++ b/src/thermo/StateModel.h @@ -120,7 +120,7 @@ class StateModel */ virtual void setState( const double* const p_mass, const double* const p_energy, - const int vars = 0, const double Tv_old = 0, const double T_old = 0) = 0; + const int vars = 0) = 0; /** * Sets the current magnitude of the magnetic field in teslas. From 58e4e40df70e87a63d8aa4b8cbc3b0a23dbf9fc7 Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Mon, 3 Feb 2025 16:00:55 -0700 Subject: [PATCH 07/23] Revert "Builds successfully. Removes T_old, Tv_old" This reverts commit 2b41465161e0570c8cba562ca36482a52911644a. --- src/thermo/Thermodynamics.cpp | 2 +- src/thermo/Thermodynamics.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index ddad37ac..32106175 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -176,7 +176,7 @@ bool Thermodynamics::speciesThermoValidAtT(const size_t i, const double T) const //============================================================================== void Thermodynamics::setState( - const double* const p_v1, const double* const p_v2, const int vars) + const double* const p_v1, const double* const p_v2, const int vars, const double Tv_old, const double T_old) { mp_state->setState(p_v1, p_v2, vars); convert(X(), mp_y); diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index 9ba1b888..aad3c9fa 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -306,7 +306,7 @@ class Thermodynamics //: public StateModelUpdateHandler * used. */ void setState( - const double* const p_v1, const double* const p_v2, const int vars = 0); + const double* const p_v1, const double* const p_v2, const int vars = 0, const double Tv_old = 0, const double T_old = 0); /** * Computes the equilibrium composition of the mixture at the given fixed From c08f303f830904e0db718d1538619f4b762d3dfa Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Mon, 3 Feb 2025 16:02:35 -0700 Subject: [PATCH 08/23] Revert "Fixed redefinition of default arguments" This reverts commit bd2f4824386866a39e7b853618e227a843e799dc. --- src/thermo/Thermodynamics.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index 32106175..5396e9e0 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -176,7 +176,7 @@ bool Thermodynamics::speciesThermoValidAtT(const size_t i, const double T) const //============================================================================== void Thermodynamics::setState( - const double* const p_v1, const double* const p_v2, const int vars, const double Tv_old, const double T_old) + const double* const p_v1, const double* const p_v2, const int vars, const double Tv_old = 0, const double T_old = 0) { mp_state->setState(p_v1, p_v2, vars); convert(X(), mp_y); From b43bbb08b34b34669c318ae1a19257865f3bfda9 Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Mon, 3 Feb 2025 16:02:37 -0700 Subject: [PATCH 09/23] Revert "Updated function headers to include Tv_old and T_old to protect against NAN in SU2" This reverts commit 382b64dd55577ab45d49070ece6409733427adf9. --- src/thermo/EquilStateModel.cpp | 2 +- src/thermo/Thermodynamics.cpp | 2 +- src/thermo/Thermodynamics.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/thermo/EquilStateModel.cpp b/src/thermo/EquilStateModel.cpp index cecda860..c8155cb9 100644 --- a/src/thermo/EquilStateModel.cpp +++ b/src/thermo/EquilStateModel.cpp @@ -272,7 +272,7 @@ class EquilTPStateModel : public EquilStateModel * Sets the state of the equilibrium mixture given temperature and pressure. */ virtual void setState( - const double* const p_T, const double* const p_P, const int vars = 1, const double Tv_old = 0, const double T_old = 0) + const double* const p_T, const double* const p_P, const int vars = 1) { //assert(vars == 1); EquilStateModel::setState(p_P, p_T, 1); diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index 5396e9e0..ddad37ac 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -176,7 +176,7 @@ bool Thermodynamics::speciesThermoValidAtT(const size_t i, const double T) const //============================================================================== void Thermodynamics::setState( - const double* const p_v1, const double* const p_v2, const int vars, const double Tv_old = 0, const double T_old = 0) + const double* const p_v1, const double* const p_v2, const int vars) { mp_state->setState(p_v1, p_v2, vars); convert(X(), mp_y); diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index aad3c9fa..9ba1b888 100644 --- a/src/thermo/Thermodynamics.h +++ b/src/thermo/Thermodynamics.h @@ -306,7 +306,7 @@ class Thermodynamics //: public StateModelUpdateHandler * used. */ void setState( - const double* const p_v1, const double* const p_v2, const int vars = 0, const double Tv_old = 0, const double T_old = 0); + const double* const p_v1, const double* const p_v2, const int vars = 0); /** * Computes the equilibrium composition of the mixture at the given fixed From e3b5c196fa7960aa8826743b8a5239d37e06ce0e Mon Sep 17 00:00:00 2001 From: Aaron Larsen Date: Mon, 3 Feb 2025 16:03:11 -0700 Subject: [PATCH 10/23] Revert "New file from Catarina that helps with NAN in SU2 coupling, See 1/20/25 email" This reverts commit 87e1c1bbe7c18ea6f3a11dce8158308e433a7489. --- src/thermo/ChemNonEqTTvStateModel.cpp | 46 +++++++++++++++++---------- 1 file changed, 29 insertions(+), 17 deletions(-) diff --git a/src/thermo/ChemNonEqTTvStateModel.cpp b/src/thermo/ChemNonEqTTvStateModel.cpp index c94e2b68..3a30acc3 100644 --- a/src/thermo/ChemNonEqTTvStateModel.cpp +++ b/src/thermo/ChemNonEqTTvStateModel.cpp @@ -74,19 +74,22 @@ class ChemNonEqTTvStateModel : public StateModel * following: * 0: {species densities}, {total energy density, vib. energy density} * 1: {species densities} and {T, Tv} + * 2: {species mass fractions} and {P, T, Tv} array */ void setState( const double* const p_mass, const double* const p_energy, - const int vars = 0, const double Tv_old = 0, const double T_old = 0) + const int vars = 0) { const int ns = m_thermo.nSpecies(); // Compute the species concentrations which are used through out this - // method regardless of variable set + // method regardless of variable set. Note that in the case where + // the p_mass vector is mass fractions, this is just dividing by the + // molecular weights as the first step in converting to mole fractions. double conc = 0.0; for (int i = 0; i < ns; ++i){ // Check that species densities are at least positive - assert(p_mass[i] >= 0.0); + assert(p_mass[i] >= 0.0 && "Species is negative!"); mp_X[i] = std::max(p_mass[i] / m_thermo.speciesMw(i), 0.0); conc += mp_X[i]; } @@ -98,38 +101,50 @@ class ChemNonEqTTvStateModel : public StateModel switch (vars) { case 0: { - solveEnergies(p_mass, p_energy, Tv_old, T_old); + solveEnergies(p_mass, p_energy); break; } case 1: // Check that temperature is at least positive - assert(p_energy[0] > 0.0); - assert(p_energy[1] > 0.0); + assert(p_energy[0] > 0.0 && "T is negative!"); + assert(p_energy[1] > 0.0 && "Tv is negative!"); m_T = p_energy[0]; m_Tv = p_energy[1]; break; + case 2: + // Check temperature and pressure are positive + assert(p_energy[0] > 0.0 && "P is negative!"); + assert(p_energy[1] > 0.0 && "T is negative!"); + assert(p_energy[2] > 0.0 && "Tv is negative!"); + m_P = p_energy[0]; + m_T = p_energy[1]; + m_Tv = p_energy[2]; + break; default: throw InvalidInputError("variable set", vars) - << "This variable-set is not implemented in ChemNonEqStateModel" + << "This variable-set is not implemented in ChemNonEqTTv" << ". Possible variable-sets are:\n" << " 0: (species densities, static energy density)\n" - << " 1: (species densities, T and Tv)"; + << " 1: (species densities, T and Tv)\n" + << " 2: (species mass fractions, P, T and Tv)"; } // Set Tr, Tel, and Te m_Tr = m_T; m_Tel = m_Te = m_Tv; - // Compute the pressure and species mole fractions from T and rho_i + // Compute the species mole fractions for (int i = 0; i < ns; ++i) mp_X[i] /= conc; - m_P = conc * RU * (m_T + (m_Tv - m_T) * elec / conc); + + // Compute the pressure from T and rho_i + if (vars != 2) + m_P = conc * RU * (m_T + (m_Tv - m_T) * elec / conc); } + /** - * Add proper - * 0: {species densities}, {total energy density, vib. energy density} - * 1: {species densities} and {T, Tv} + * Initialize transfer models */ void initializeTransferModel(Mutation::Mixture& mix) { @@ -239,7 +254,7 @@ class ChemNonEqTTvStateModel : public StateModel * @param p_rhoe - total and internal mass energy */ - void solveEnergies(const double* const p_rhoi, const double* const p_rhoe, const double Tv_old = 0, const double T_old = 0) + void solveEnergies(const double* const p_rhoi, const double* const p_rhoe) { // This line is required to have perfect numerical jacobians... //m_T = m_Tv = 500.0; @@ -280,9 +295,6 @@ class ChemNonEqTTvStateModel : public StateModel f = e - emix; } - if (m_T != m_T) m_T = T_old; - if (m_Tv != m_Tv) m_Tv = Tv_old; - if (i == imax) std::cout << "Warning, didn't converge temperatures: f = " << f.norm() << std::endl; } From 405b6d4da25a1540b968bacdebf7422d7b58eb86 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Wed, 9 Apr 2025 14:16:49 -0400 Subject: [PATCH 11/23] Added some variables and functions to caluclate the temperature-dependent chemistry jacobians. --- src/kinetics/JacobianManager.cpp | 34 +++++++ src/kinetics/JacobianManager.h | 17 ++++ src/kinetics/Kinetics.cpp | 69 ++++++++++++- src/kinetics/Kinetics.h | 13 +++ src/kinetics/RateLawGroup.h | 164 ++++++++++++++++++++++++++++++- src/kinetics/RateLaws.h | 50 ++++++++++ src/kinetics/RateManager.cpp | 41 +++++++- src/kinetics/RateManager.h | 21 ++++ 8 files changed, 401 insertions(+), 8 deletions(-) diff --git a/src/kinetics/JacobianManager.cpp b/src/kinetics/JacobianManager.cpp index 4f718a0e..683fcf4e 100644 --- a/src/kinetics/JacobianManager.cpp +++ b/src/kinetics/JacobianManager.cpp @@ -315,6 +315,40 @@ void JacobianManager::computeJacobian( } } +//============================================================================== + +void JacobianManager::computeJacobianT( + const double* const p_dropf, const double* const p_dropb, + double* const Tjac) const +{ + const size_t ns = m_thermo.nSpecies(); + const size_t nr = m_reactions.size(); + const size_t nt = m_thermo.nEnergyEqns(); + + double* m_t; + double* p_t = new double[nt]; + + std::fill(p_t, p_t + nt, 0.); + + m_thermo.getTemperatures(p_t); + + // Make sure we are starting with a clean slate + std::fill(Tjac, Tjac+ns*nt, 0.0); + + +} +/* +//============================================================================== + +void JacobianManager::computeTReactionDerivative( + const size_t nr, const size_t nt, const double* const p_t, + double* const p_dTfdT, double*const p_dTbdT) const +{ + //for(size_t i = 0; i < nr; ++i) { + + +} +*/ //============================================================================== } // namespace Kinetics diff --git a/src/kinetics/JacobianManager.h b/src/kinetics/JacobianManager.h index 1f8597a7..0d2a65d0 100644 --- a/src/kinetics/JacobianManager.h +++ b/src/kinetics/JacobianManager.h @@ -31,6 +31,7 @@ #include "Thermodynamics.h" #include +#include "RateLaws.h" // Added by RSCD namespace Mutation { namespace Kinetics { @@ -517,7 +518,23 @@ class JacobianManager void computeJacobian( const double* const kf, const double* const kb, const double* const conc, double* const sjac) const; + + /** + * Computes the temperature-dependent source Jacobian elements. + * Work under progress - RSCD + */ + void computeJacobianT( + const double* const p_ropf, const double* const p_ropb, + double* const Tjac) const; + /** + * Computes the reaction temperature derivatives. + * Work under progress - RSCD + */ +/* void computeTReactionDerivative( + const size_t nr, const size_t nt, const double* const p_t, + double* const p_dTfdT, double* const p_dTbdT) const; +*/ private: /** diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 9cdd55a1..9207e994 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -54,7 +54,9 @@ Kinetics::Kinetics( mp_ropf(NULL), mp_ropb(NULL), mp_rop(NULL), - mp_wdot(NULL) + mp_wdot(NULL), + mp_dropf(NULL), + mp_dropb(NULL) { if (mechanism == "none") return; @@ -104,6 +106,10 @@ 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; } //============================================================================== @@ -162,7 +168,8 @@ 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()]; } //============================================================================== @@ -279,6 +286,45 @@ void Kinetics::backwardRatesOfProgress( m_thirdbodies.multiplyThirdbodies(p_conc, p_ropb); } +//============================================================================== + +void Kinetics::forwardRateOfProgressDerivatives(double* const mp_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) + mp_dropf[i] = mp_ropf[i]*dkfdT[i]; + +} + + +//============================================================================== + +void Kinetics::backwardRateOfProgressDerivatives(double* const mp_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) + mp_dropb[i] = mp_ropb[i]*dkbdT[i]; + + +} + + /* //============================================================================== @@ -387,6 +433,25 @@ void Kinetics::jacobianRho(double* const p_jac) m_jacobian.computeJacobian(mp_ropf, mp_ropb, mp_rop, p_jac); } +//============================================================================== + +void Kinetics::jacobianT(double* const p_jacT, double* const p_jacTv) +{ + // Special case of no reactions + if (nReactions() == 0) { + std::fill(p_jacT, p_jacT + m_thermo.nSpecies(), 0); + std::fill(p_jacTv, p_jacTv + m_thermo.nSpecies(), 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, p_jacT); +} + //============================================================================== } // namespace Kinetics diff --git a/src/kinetics/Kinetics.h b/src/kinetics/Kinetics.h index eb100139..2118acab 100644 --- a/src/kinetics/Kinetics.h +++ b/src/kinetics/Kinetics.h @@ -194,6 +194,17 @@ class Kinetics */ void jacobianRho(double* const p_jac); + /** + * Fills the matrix p_jacT with the temperature jacobians for chemistry source term + * \f[ + * J_{ij} = \frac{\partial \dot{\omega}_i}{\partial T_j} + * \f] + * Work in progress and more details will be added soon - RSCD + */ + void jacobianT(double* const p_jacT, double* const p_jacTv); + + void forwardRateOfProgressDerivatives(double* const p_dropf); + void backwardRateOfProgressDerivatives(double* const p_dropb); /** * Returns the change in some species quantity across each reaction. */ @@ -237,6 +248,8 @@ class Kinetics double* mp_ropb; double* mp_rop; double* mp_wdot; + double* mp_dropf; + double* mp_dropb; }; diff --git a/src/kinetics/RateLawGroup.h b/src/kinetics/RateLawGroup.h index 78b7b679..83c48e48 100644 --- a/src/kinetics/RateLawGroup.h +++ b/src/kinetics/RateLawGroup.h @@ -93,6 +93,18 @@ 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; + + virtual void dTreacdTtr( + const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdT) = 0; + + virtual void dTreacdTve( + const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdTv) = 0; + /** * Computes \Delta G / RT for this rate law group and subtracts these values * for each of the reactions in this group. @@ -103,12 +115,21 @@ class RateLawGroup const double val = std::log(ONEATM / (RU * m_t)); for (int 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(size_t ns, double* const p_dKeqdT, double* const p_dkdT) const + { + m_reacs.decrReactions(p_dKeqdT, p_dkdT); + m_prods.incrReactions(p_dKeqdT, p_dkdT); + } protected: @@ -116,6 +137,7 @@ class RateLawGroup /// in the lnk() function) double m_t; double m_last_t; + const char* m_name; /// Stores the reactants for reactions that will use this rate law for the /// reverse direction @@ -172,6 +194,83 @@ 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 lnT = std::log(m_t); + 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; + } + + /** + * Evaluates all of the reaction temperature derivatives in the group + */ + virtual void dTreacdTtr( + const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdT) + { + // Determine the reaction temperature for this group + m_t = TSelectorType().getT(thermo.state()); + + // Determine the TSelector type for this group + m_name = TSelectorType().getName(thermo.state()); + +/* //Number of energy equations + const size_t nt = thermo.nEnergyEqns(); + + //Temperatures of the state + double* p_t = new double[nt](); + //std::fill(p_t, p_t + nt, 0.); + thermo.getTemperatures(p_t); +*/ + for (size_t i = 0; i < m_rates.size(); ++i) { + const std::pair& rate = m_rates[i]; + p_dTfdT[rate.first] = rate.second.dTrxndTtr(m_name, m_t, nt, p_t); + } + + } + + /** + * Evaluates all of the reaction temperature derivatives in the group + */ + virtual void dTreacdTve( + const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdTv) + { + // Determine the reaction temperature for this group + m_t = TSelectorType().getT(thermo.state()); + + // Determine the TSelector type for this group + m_name = TSelectorType().getName(thermo.state()); + +/* //Number of energy equations + const size_t nt = thermo.nEnergyEqns(); + + //Temperatures of the state + double* p_t = new double[nt](); + //std::fill(p_t, p_t + nt, 0.); + thermo.getTemperatures(p_t); +*/ + for (size_t i = 0; i < m_rates.size(); ++i) { + const std::pair& rate = m_rates[i]; + p_dTfdTv[rate.first] = rate.second.dTrxndTve(m_name, m_t, nt, p_t); + } + } + private: /// vector of rates to evaluate @@ -259,6 +358,19 @@ class RateLawGroupCollection iter->second->lnk(p_state, p_lnk); } + /** + * Computes the rate coefficients 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 +387,56 @@ class RateLawGroupCollection } } + /** + * Computes the temperature derivative of keq as follows: + * \f[ + * \frac{1}{keq} \frac{\partial keq}{\partial T} = + * \frac{1}{T} \sum_{i=1}^{ns} {\delta\nu_{ij}} \left(\frac{h_i}{Ru T}-1\right) + * \f] + * More details will be added soon - RSCD + */ + 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.speciesHOverRT(p_group->getT(), p_dKeqdT); + for(size_t i = 0; i < ns; ++i) { + p_dKeqdT[i] = 1. - p_dKeqdT[i]; + p_dKeqdT[i] /= (p_group->getT()); + } + p_group->derivativeKeq(ns, p_dKeqdT, p_dkdT); + } + } + + /** + * Computes the reaction temperature derivatives with respect Ttr and Tve + * More details will be added soon - RSCD + */ + void derivativeTrxnTtr( + const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdT) + { + GroupMap::iterator iter = m_group_map.begin(); + for ( ; iter != m_group_map.end(); ++iter) + iter->second->dTreacdTtr(thermo, nt, p_t, p_dTfdT); + + } + + /** + * Computes the reaction temperature derivatives with respect Ttr and Tve + * More details will be added soon - RSCD + */ + void derivativeTrxnTve( + const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdTv) + { + GroupMap::iterator iter = m_group_map.begin(); + for ( ; iter != m_group_map.end(); ++iter) + iter->second->dTreacdTve(thermo, nt, p_t, p_dTfdTv); + + } + private: /// Collection of RateLawGroup objects diff --git a/src/kinetics/RateLaws.h b/src/kinetics/RateLaws.h index 8f9ac6b8..4321cfc0 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,6 +93,52 @@ class Arrhenius : public RateLaw double T() const { return m_temp; } + + inline double dTrxndTtr(const char* m_name, const double m_t, const size_t nt, const double* const p_t) const { +/* std::string selector(m_name); + if (selector == "TSelector") { + return (m_t/p_t[0]); + } + else if (selector == "TeSelector") { + return 0.; + } + else if (selector == "ParkSelector") { + return (0.5*m_t/p_t[0]); + } + else { + throw InvalidInputError("TSelector type", selector) + << "TSelector type mentioned is not implemented."; + } +*/ + std::string selector(m_name); + if (m_t == p_t[0]) return 1.; + else if (m_t == p_t[1]) return 0.; + else return (0.5*m_t/p_t[0]); + } + + inline double dTrxndTve(const char* m_name, const double m_t, const size_t nt, const double* const p_t) const { + switch (nt) { + case 1: + return 0.0; + case 2: + if (std::string(m_name) == "ParkSelector") { + return (0.5*m_t/p_t[1]); + } + else if (std::string(m_name) == "TSelector") { + return 0.0; + } + else if (std::string(m_name) == "TeSelector") { + return 1.0; + } + else { + throw InvalidInputError("TSelector type", m_name) + << "TSelector type mentioned is not implemented."; + } + break; + default: + return 0.0; + } + } private: diff --git a/src/kinetics/RateManager.cpp b/src/kinetics/RateManager.cpp index 40d4460b..37985e4b 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,29 @@ 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), + mp_dTfdT(NULL), mp_dTfdTv(NULL), mp_dTbdT(NULL), mp_dTbdTv(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) + for (size_t i = 0; i < m_nr; ++i) { addReaction(i, reactions[i]); - + //std::cout<<"reaction:"<<'\t'<(rxn, reaction); @@ -210,7 +223,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 +238,34 @@ 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 rate coefficients temperature derivatives + 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 temperature derivatives of equilibrium constants + m_rate_groups.derivativeKeq(thermo, mp_dKeqdT, mp_dkbdT); + //Number of energy equations + const size_t nt = thermo.nEnergyEqns(); + + //Temperatures of the state + double* p_t = new double[nt]; + std::fill(p_t, p_t + nt, 0.); + thermo.getTemperatures(p_t); + + // Evaluate the reaction temperature derivatives + m_rate_groups.derivativeTrxnTtr(thermo, nt, p_t, mp_dTfdT); + m_rate_groups.derivativeTrxnTve(thermo, nt, p_t, mp_dTfdTv); } //============================================================================== diff --git a/src/kinetics/RateManager.h b/src/kinetics/RateManager.h index 5d93148c..44624e56 100644 --- a/src/kinetics/RateManager.h +++ b/src/kinetics/RateManager.h @@ -79,6 +79,12 @@ class RateManager return m_irr; } + /* need to add once it works - RSCD + */ + + const double* const dkfdT() {return mp_dkfdT; } + const double* const dkbdT() {return mp_dkbdT; } + private: /** @@ -122,6 +128,21 @@ 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; + + /// Storage for reaction temperature derivatives + double* mp_dTfdT; + double* mp_dTfdTv; + double* mp_dTbdT; + double* mp_dTbdTv; + /// Stores the indices for which the forward and reverse temperature /// evaluations are equal std::vector m_to_copy; From c307612ef5afb5679231d4013530ec11b492d793 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Tue, 15 Apr 2025 21:25:10 -0400 Subject: [PATCH 12/23] Finished adding necessary functions and variables that will help in calculating temperature-related chemistry jacobians. --- src/kinetics/JacobianManager.cpp | 190 +++++++++++++++++++++++++--- src/kinetics/JacobianManager.h | 35 +++-- src/kinetics/Kinetics.cpp | 64 ++++++++-- src/kinetics/Kinetics.h | 43 ++++++- src/kinetics/RateLawGroup.h | 107 ++-------------- src/kinetics/RateLaws.h | 46 ------- src/kinetics/RateManager.cpp | 23 +--- src/kinetics/RateManager.h | 26 ++-- src/kinetics/StoichiometryManager.h | 2 +- 9 files changed, 314 insertions(+), 222 deletions(-) diff --git a/src/kinetics/JacobianManager.cpp b/src/kinetics/JacobianManager.cpp index 683fcf4e..7bdcb3d5 100644 --- a/src/kinetics/JacobianManager.cpp +++ b/src/kinetics/JacobianManager.cpp @@ -293,6 +293,156 @@ 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::computeJacobian( const double* const kf, const double* const kb, const double* const conc, double* const sjac) const @@ -319,36 +469,42 @@ void JacobianManager::computeJacobian( void JacobianManager::computeJacobianT( const double* const p_dropf, const double* const p_dropb, - double* const Tjac) const + const std::vector& reactions, double* const p_dropT) const { - const size_t ns = m_thermo.nSpecies(); - const size_t nr = m_reactions.size(); const size_t nt = m_thermo.nEnergyEqns(); - - double* m_t; double* p_t = new double[nt]; - + std::vector > dTrxndT; + std::fill(p_t, p_t + nt, 0.); - m_thermo.getTemperatures(p_t); + + computeTReactionTDerivatives(reactions, p_t, dTrxndT); - // Make sure we are starting with a clean slate - std::fill(Tjac, Tjac+ns*nt, 0.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::computeTReactionDerivative( - const size_t nr, const size_t nt, const double* const p_t, - double* const p_dTfdT, double*const p_dTbdT) const +void JacobianManager::computeJacobianTv( + const double* const p_dropf, const double* const p_dropb, + const std::vector& reactions, double* const p_dropTv) const { - //for(size_t i = 0; i < nr; ++i) { - + const size_t nt = m_thermo.nEnergyEqns(); + double* p_t = new double[nt]; + std::vector > dTrxndTv; + + std::fill(p_t, p_t + nt, 0.); + m_thermo.getTemperatures(p_t); + + computeTReactionTvDerivatives(reactions, p_t, dTrxndTv); + + 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 0d2a65d0..247a478f 100644 --- a/src/kinetics/JacobianManager.h +++ b/src/kinetics/JacobianManager.h @@ -31,7 +31,6 @@ #include "Thermodynamics.h" #include -#include "RateLaws.h" // Added by RSCD namespace Mutation { namespace Kinetics { @@ -520,21 +519,37 @@ class JacobianManager const double* const conc, double* const sjac) const; /** - * Computes the temperature-dependent source Jacobian elements. - * Work under progress - RSCD + * Computes the translational-temperature-dependent + * chemistry Jacobian elements. */ void computeJacobianT( const double* const p_ropf, const double* const p_ropb, - double* const Tjac) const; + const std::vector& reactions, double* const p_dropT) const; /** - * Computes the reaction temperature derivatives. - * Work under progress - RSCD + * Computes the vibrational-temperature-dependent + * chemistry Jacobian elements. */ -/* void computeTReactionDerivative( - const size_t nr, const size_t nt, const double* const p_t, - double* const p_dTfdT, double* const p_dTbdT) const; -*/ + 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: /** diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 9207e994..8b014248 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -56,7 +56,9 @@ Kinetics::Kinetics( mp_rop(NULL), mp_wdot(NULL), mp_dropf(NULL), - mp_dropb(NULL) + mp_dropb(NULL), + mp_dropT(NULL), + mp_dropTv(NULL) { if (mechanism == "none") return; @@ -110,6 +112,10 @@ Kinetics::~Kinetics() delete [] mp_dropf; if (mp_dropb != NULL) delete [] mp_dropb; + if (mp_dropT != NULL) + delete [] mp_dropT; + if (mp_dropTv != NULL) + delete [] mp_dropTv; } //============================================================================== @@ -170,6 +176,8 @@ void Kinetics::closeReactions(const bool validate_mechanism) mp_wdot = new double [m_thermo.nSpecies()]; mp_dropf = new double [nReactions()]; mp_dropb = new double [nReactions()]; + mp_dropT = new double [nReactions()]; + mp_dropTv = new double [nReactions()]; } //============================================================================== @@ -299,12 +307,10 @@ void Kinetics::forwardRateOfProgressDerivatives(double* const mp_dropf) forwardRatesOfProgress(mp_ropf); - for (int i = 0; i < nReactions(); ++i) + for (int i = 0; i < nReactions(); ++i) mp_dropf[i] = mp_ropf[i]*dkfdT[i]; - } - //============================================================================== void Kinetics::backwardRateOfProgressDerivatives(double* const mp_dropb) @@ -315,13 +321,11 @@ void Kinetics::backwardRateOfProgressDerivatives(double* const mp_dropb) mp_rates->update(m_thermo); ArrayXd dkbdT = Map(mp_rates->dkbdT(), nReactions()); - + backwardRatesOfProgress(mp_ropb); for(int i=0; i < nReactions(); ++i) mp_dropb[i] = mp_ropb[i]*dkbdT[i]; - - } @@ -435,21 +439,53 @@ void Kinetics::jacobianRho(double* const p_jac) //============================================================================== -void Kinetics::jacobianT(double* const p_jacT, double* const p_jacTv) +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) { - std::fill(p_jacT, p_jacT + m_thermo.nSpecies(), 0); - std::fill(p_jacTv, p_jacTv + m_thermo.nSpecies(), 0); - return; - } + 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.computeJacobianT(mp_dropf, mp_dropb, p_jacT); + 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 2118acab..68bc1ab6 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 @@ -197,14 +217,25 @@ class Kinetics /** * Fills the matrix p_jacT with the temperature jacobians for chemistry source term * \f[ - * J_{ij} = \frac{\partial \dot{\omega}_i}{\partial T_j} + * J_{i} = \frac{\partial \dot{\omega}_i}{\partial T_{tr}} * \f] - * Work in progress and more details will be added soon - RSCD + * 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, double* const p_jacTv); + 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); - void forwardRateOfProgressDerivatives(double* const p_dropf); - void backwardRateOfProgressDerivatives(double* const p_dropb); /** * Returns the change in some species quantity across each reaction. */ @@ -250,6 +281,8 @@ class Kinetics double* mp_wdot; double* mp_dropf; double* mp_dropb; + double* mp_dropT; + double* mp_dropTv; }; diff --git a/src/kinetics/RateLawGroup.h b/src/kinetics/RateLawGroup.h index 83c48e48..42c839af 100644 --- a/src/kinetics/RateLawGroup.h +++ b/src/kinetics/RateLawGroup.h @@ -99,12 +99,6 @@ class RateLawGroup virtual void invkdkdT( const Thermodynamics::StateModel* const p_state, double* const p_dkdT) = 0; - virtual void dTreacdTtr( - const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdT) = 0; - - virtual void dTreacdTve( - const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdTv) = 0; - /** * Computes \Delta G / RT for this rate law group and subtracts these values * for each of the reactions in this group. @@ -113,7 +107,7 @@ 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 @@ -125,7 +119,7 @@ class RateLawGroup * 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(size_t ns, double* const p_dKeqdT, double* const p_dkdT) const + 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); @@ -137,7 +131,6 @@ class RateLawGroup /// in the lnk() function) double m_t; double m_last_t; - const char* m_name; /// Stores the reactants for reactions that will use this rate law for the /// reverse direction @@ -205,7 +198,6 @@ class RateLawGroup1T : public RateLawGroup // Update only if the temperature has changed //if (std::abs(m_t - m_last_t) > 1.0e-10) { - const double lnT = std::log(m_t); const double invT = 1.0 / m_t; for (size_t i = 0; i < m_rates.size(); ++i) { @@ -218,59 +210,6 @@ class RateLawGroup1T : public RateLawGroup m_last_t = m_t; } - /** - * Evaluates all of the reaction temperature derivatives in the group - */ - virtual void dTreacdTtr( - const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdT) - { - // Determine the reaction temperature for this group - m_t = TSelectorType().getT(thermo.state()); - - // Determine the TSelector type for this group - m_name = TSelectorType().getName(thermo.state()); - -/* //Number of energy equations - const size_t nt = thermo.nEnergyEqns(); - - //Temperatures of the state - double* p_t = new double[nt](); - //std::fill(p_t, p_t + nt, 0.); - thermo.getTemperatures(p_t); -*/ - for (size_t i = 0; i < m_rates.size(); ++i) { - const std::pair& rate = m_rates[i]; - p_dTfdT[rate.first] = rate.second.dTrxndTtr(m_name, m_t, nt, p_t); - } - - } - - /** - * Evaluates all of the reaction temperature derivatives in the group - */ - virtual void dTreacdTve( - const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdTv) - { - // Determine the reaction temperature for this group - m_t = TSelectorType().getT(thermo.state()); - - // Determine the TSelector type for this group - m_name = TSelectorType().getName(thermo.state()); - -/* //Number of energy equations - const size_t nt = thermo.nEnergyEqns(); - - //Temperatures of the state - double* p_t = new double[nt](); - //std::fill(p_t, p_t + nt, 0.); - thermo.getTemperatures(p_t); -*/ - for (size_t i = 0; i < m_rates.size(); ++i) { - const std::pair& rate = m_rates[i]; - p_dTfdTv[rate.first] = rate.second.dTrxndTve(m_name, m_t, nt, p_t); - } - } - private: /// vector of rates to evaluate @@ -359,8 +298,10 @@ class RateLawGroupCollection } /** - * Computes the rate coefficients in this collection and stores them in - * the vector at the index corresponding to their respective reaction. + * 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) @@ -389,11 +330,9 @@ class RateLawGroupCollection /** * Computes the temperature derivative of keq as follows: - * \f[ - * \frac{1}{keq} \frac{\partial keq}{\partial T} = - * \frac{1}{T} \sum_{i=1}^{ns} {\delta\nu_{ij}} \left(\frac{h_i}{Ru T}-1\right) - * \f] - * More details will be added soon - RSCD + * \f$ \frac{1}{keq} \frac{\partial keq}{\partial T} = + * \frac{1}{T} \sum_{i=1}^{ns} {\delta\nu_{ij}} \left(\frac{h_i}{Ru T}-1\right) + * \f$ */ void derivativeKeq( const Thermodynamics::Thermodynamics& thermo, double* const p_dKeqdT, double* const p_dkdT) @@ -407,35 +346,9 @@ class RateLawGroupCollection p_dKeqdT[i] = 1. - p_dKeqdT[i]; p_dKeqdT[i] /= (p_group->getT()); } - p_group->derivativeKeq(ns, p_dKeqdT, p_dkdT); + p_group->derivativeKeq(p_dKeqdT, p_dkdT); } } - - /** - * Computes the reaction temperature derivatives with respect Ttr and Tve - * More details will be added soon - RSCD - */ - void derivativeTrxnTtr( - const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdT) - { - GroupMap::iterator iter = m_group_map.begin(); - for ( ; iter != m_group_map.end(); ++iter) - iter->second->dTreacdTtr(thermo, nt, p_t, p_dTfdT); - - } - - /** - * Computes the reaction temperature derivatives with respect Ttr and Tve - * More details will be added soon - RSCD - */ - void derivativeTrxnTve( - const Thermodynamics::Thermodynamics& thermo, const size_t nt, const double* const p_t, double* const p_dTfdTv) - { - GroupMap::iterator iter = m_group_map.begin(); - for ( ; iter != m_group_map.end(); ++iter) - iter->second->dTreacdTve(thermo, nt, p_t, p_dTfdTv); - - } private: diff --git a/src/kinetics/RateLaws.h b/src/kinetics/RateLaws.h index 4321cfc0..85a45458 100644 --- a/src/kinetics/RateLaws.h +++ b/src/kinetics/RateLaws.h @@ -94,52 +94,6 @@ class Arrhenius : public RateLaw return m_temp; } - inline double dTrxndTtr(const char* m_name, const double m_t, const size_t nt, const double* const p_t) const { -/* std::string selector(m_name); - if (selector == "TSelector") { - return (m_t/p_t[0]); - } - else if (selector == "TeSelector") { - return 0.; - } - else if (selector == "ParkSelector") { - return (0.5*m_t/p_t[0]); - } - else { - throw InvalidInputError("TSelector type", selector) - << "TSelector type mentioned is not implemented."; - } -*/ - std::string selector(m_name); - if (m_t == p_t[0]) return 1.; - else if (m_t == p_t[1]) return 0.; - else return (0.5*m_t/p_t[0]); - } - - inline double dTrxndTve(const char* m_name, const double m_t, const size_t nt, const double* const p_t) const { - switch (nt) { - case 1: - return 0.0; - case 2: - if (std::string(m_name) == "ParkSelector") { - return (0.5*m_t/p_t[1]); - } - else if (std::string(m_name) == "TSelector") { - return 0.0; - } - else if (std::string(m_name) == "TeSelector") { - return 1.0; - } - else { - throw InvalidInputError("TSelector type", m_name) - << "TSelector type mentioned is not implemented."; - } - break; - default: - return 0.0; - } - } - private: static std::vector sm_aunits; diff --git a/src/kinetics/RateManager.cpp b/src/kinetics/RateManager.cpp index 37985e4b..557a541e 100644 --- a/src/kinetics/RateManager.cpp +++ b/src/kinetics/RateManager.cpp @@ -117,14 +117,12 @@ 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_dkfdT(NULL), mp_dkbdT(NULL), mp_dKeqdT(NULL), - mp_dTfdT(NULL), mp_dTfdTv(NULL), mp_dTbdT(NULL), mp_dTbdTv(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]); - //std::cout<<"reaction:"<<'\t'<& reactions) mp_dkfdT = mp_gibbs + block_size; mp_dkbdT = mp_dkfdT + m_nr; mp_dKeqdT = mp_dkbdT + m_nr; - mp_dTfdT = mp_dKeqdT + block_size; - mp_dTfdTv = mp_dTfdT + block_size; - mp_dTbdT = mp_dTfdTv + m_nr; - mp_dTbdTv = mp_dTbdT + m_nr; // Initialize the arrays to zero std::fill(mp_lnkf, mp_lnkf+block_size, 0.0); @@ -238,7 +232,8 @@ 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 rate coefficients temperature derivatives + // 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 @@ -253,19 +248,9 @@ void RateManager::update(const Thermodynamics::Thermodynamics& thermo) // Subtract lnkeq(Tb) rate constants from the lnkf(Tb) to get lnkb(Tb) m_rate_groups.subtractLnKeq(thermo, mp_gibbs, mp_lnkb); - // Evaluate the temperature derivatives of equilibrium constants + // Evaluate the reaction temperature derivatives of equilibrium constants m_rate_groups.derivativeKeq(thermo, mp_dKeqdT, mp_dkbdT); - //Number of energy equations - const size_t nt = thermo.nEnergyEqns(); - //Temperatures of the state - double* p_t = new double[nt]; - std::fill(p_t, p_t + nt, 0.); - thermo.getTemperatures(p_t); - - // Evaluate the reaction temperature derivatives - m_rate_groups.derivativeTrxnTtr(thermo, nt, p_t, mp_dTfdT); - m_rate_groups.derivativeTrxnTve(thermo, nt, p_t, mp_dTfdTv); } //============================================================================== diff --git a/src/kinetics/RateManager.h b/src/kinetics/RateManager.h index 44624e56..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,11 +79,17 @@ class RateManager return m_irr; } - /* need to add once it works - RSCD - */ - - const double* const dkfdT() {return mp_dkfdT; } - const double* const dkbdT() {return mp_dkbdT; } + /** + * 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: @@ -137,12 +143,6 @@ class RateManager /// Storage for equilibrium constant temperature derivatives double* mp_dKeqdT; - /// Storage for reaction temperature derivatives - double* mp_dTfdT; - double* mp_dTfdTv; - double* mp_dTbdT; - double* mp_dTbdTv; - /// 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]]; } From dac2a5dd484d940bf22ff1d948913d847809e86d Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Fri, 6 Jun 2025 00:36:21 -0700 Subject: [PATCH 13/23] Added functions and variables that help in computing the vibrational source term jacobians and also made changes to earlier version to minimize the memory leak. --- compile.sh | 4 + data/mechanisms/air5p_Park.xml | 29 ++++++++ data/transport/collisions.xml | 2 +- src/general/Mixture.h | 16 ++++ src/kinetics/JacobianManager.cpp | 96 +++++++++++++++++++++--- src/kinetics/JacobianManager.h | 38 ++++++++++ src/kinetics/Kinetics.cpp | 71 +++++++++++++++++- src/kinetics/Kinetics.h | 37 ++++++++++ src/kinetics/RateLawGroup.h | 12 +-- src/kinetics/RateManager.cpp | 4 +- src/numerics/Functors.h | 13 ++++ src/thermo/.StateModel.h.swp | Bin 0 -> 16384 bytes src/thermo/HarmonicOscillator.cpp | 11 ++- src/thermo/HarmonicOscillator.h | 5 +- src/thermo/MultiPhaseEquilSolver.h | 2 +- src/thermo/Nasa7Polynomial.cpp | 23 ++++++ src/thermo/Nasa7Polynomial.h | 12 ++- src/thermo/Nasa9Polynomial.cpp | 25 ++++++- src/thermo/Nasa9Polynomial.h | 12 ++- src/thermo/NasaDB.h | 23 ++++++ src/thermo/RrhoDB.cpp | 100 +++++++++++++++++++++++-- src/thermo/StateModel.h | 46 +++++++++++- src/thermo/ThermoDB.h | 22 ++++++ src/thermo/Thermodynamics.cpp | 7 ++ src/thermo/Thermodynamics.h | 7 ++ src/transfer/OmegaCE.cpp | 31 ++++++++ src/transfer/OmegaCElec.cpp | 54 ++++++++++++++ src/transfer/OmegaCV.cpp | 95 +++++++++++++++++++++++- src/transfer/OmegaET.cpp | 62 +++++++++++++++- src/transfer/OmegaI.cpp | 113 ++++++++++++++++++++++++++++- src/transfer/OmegaVT.cpp | 47 ++++++++++++ src/transfer/TransferModel.h | 4 + 32 files changed, 977 insertions(+), 46 deletions(-) create mode 100755 compile.sh create mode 100644 data/mechanisms/air5p_Park.xml create mode 100644 src/thermo/.StateModel.h.swp diff --git a/compile.sh b/compile.sh new file mode 100755 index 00000000..9a1bbad6 --- /dev/null +++ b/compile.sh @@ -0,0 +1,4 @@ +#!/bin/bash +rm -r build; mkdir build; cd build +cmake -DCMAKE_INSTALL_PREFIX=/Users/raghava/Desktop/CHANL_HPC/Mutationpp/install .. +make -j 20 install diff --git a/data/mechanisms/air5p_Park.xml b/data/mechanisms/air5p_Park.xml new file mode 100644 index 00000000..584e7a84 --- /dev/null +++ b/data/mechanisms/air5p_Park.xml @@ -0,0 +1,29 @@ + + + + + + + + + N2:0.2333, NO:0.2333, O2:0.2333, N:1.0, O:1.0 + + + + + + + + + + + + + + + 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 7bdcb3d5..aa770925 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, @@ -443,6 +494,31 @@ void JacobianManager::computeTReactionTvDerivatives( //============================================================================== +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 @@ -454,15 +530,15 @@ 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); - } } //============================================================================== @@ -472,13 +548,12 @@ void JacobianManager::computeJacobianT( const std::vector& reactions, double* const p_dropT) const { const size_t nt = m_thermo.nEnergyEqns(); - double* p_t = new double[nt]; std::vector > dTrxndT; - std::fill(p_t, p_t + nt, 0.); - m_thermo.getTemperatures(p_t); + std::fill(mp_T, mp_T + nt, 0.); + m_thermo.getTemperatures(mp_T); - computeTReactionTDerivatives(reactions, p_t, dTrxndT); + computeTReactionTDerivatives(reactions, mp_T, dTrxndT); for (int i = 0; i < reactions.size(); ++i) p_dropT[i] = p_dropf[i]*dTrxndT[i].first - p_dropb[i]*dTrxndT[i].second; @@ -492,13 +567,12 @@ void JacobianManager::computeJacobianTv( const std::vector& reactions, double* const p_dropTv) const { const size_t nt = m_thermo.nEnergyEqns(); - double* p_t = new double[nt]; std::vector > dTrxndTv; - std::fill(p_t, p_t + nt, 0.); - m_thermo.getTemperatures(p_t); + std::fill(mp_T, mp_T + nt, 0.); + m_thermo.getTemperatures(mp_T); - computeTReactionTvDerivatives(reactions, p_t, dTrxndTv); + computeTReactionTvDerivatives(reactions, mp_T, dTrxndTv); for (int i = 0; i < reactions.size(); ++i) p_dropTv[i] = p_dropf[i]*dTrxndTv[i].first - p_dropb[i]*dTrxndTv[i].second; diff --git a/src/kinetics/JacobianManager.h b/src/kinetics/JacobianManager.h index 247a478f..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) @@ -518,6 +547,13 @@ class JacobianManager 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. @@ -576,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 8b014248..04b46926 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -57,6 +57,7 @@ Kinetics::Kinetics( mp_wdot(NULL), mp_dropf(NULL), mp_dropb(NULL), + mp_dropRho(NULL), mp_dropT(NULL), mp_dropTv(NULL) { @@ -112,6 +113,8 @@ Kinetics::~Kinetics() 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) @@ -176,6 +179,7 @@ void Kinetics::closeReactions(const bool validate_mechanism) 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()]; } @@ -296,7 +300,7 @@ void Kinetics::backwardRatesOfProgress( //============================================================================== -void Kinetics::forwardRateOfProgressDerivatives(double* const mp_dropf) +void Kinetics::forwardRateOfProgressDerivatives(double* const p_dropf) { if(nReactions() == 0) return; @@ -308,12 +312,12 @@ void Kinetics::forwardRateOfProgressDerivatives(double* const mp_dropf) forwardRatesOfProgress(mp_ropf); for (int i = 0; i < nReactions(); ++i) - mp_dropf[i] = mp_ropf[i]*dkfdT[i]; + p_dropf[i] = mp_ropf[i]*dkfdT[i]; } //============================================================================== -void Kinetics::backwardRateOfProgressDerivatives(double* const mp_dropb) +void Kinetics::backwardRateOfProgressDerivatives(double* const p_dropb) { if (nReactions() == 0) return; @@ -325,7 +329,7 @@ void Kinetics::backwardRateOfProgressDerivatives(double* const mp_dropb) backwardRatesOfProgress(mp_ropb); for(int i=0; i < nReactions(); ++i) - mp_dropb[i] = mp_ropb[i]*dkbdT[i]; + p_dropb[i] = mp_ropb[i]*dkbdT[i]; } @@ -417,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 diff --git a/src/kinetics/Kinetics.h b/src/kinetics/Kinetics.h index 68bc1ab6..7c07c0f0 100644 --- a/src/kinetics/Kinetics.h +++ b/src/kinetics/Kinetics.h @@ -201,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[ @@ -281,6 +317,7 @@ class Kinetics 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 42c839af..33d51c3c 100644 --- a/src/kinetics/RateLawGroup.h +++ b/src/kinetics/RateLawGroup.h @@ -330,9 +330,7 @@ class RateLawGroupCollection /** * Computes the temperature derivative of keq as follows: - * \f$ \frac{1}{keq} \frac{\partial keq}{\partial T} = - * \frac{1}{T} \sum_{i=1}^{ns} {\delta\nu_{ij}} \left(\frac{h_i}{Ru T}-1\right) - * \f$ + * \f$ \frac{1}{keq} \frac{\partial keq}{\partial T} \f$ */ void derivativeKeq( const Thermodynamics::Thermodynamics& thermo, double* const p_dKeqdT, double* const p_dkdT) @@ -341,11 +339,9 @@ class RateLawGroupCollection GroupMap::iterator iter = m_group_map.begin(); for ( ; iter != m_group_map.end(); ++iter) { const RateLawGroup* p_group = iter->second; - thermo.speciesHOverRT(p_group->getT(), p_dKeqdT); - for(size_t i = 0; i < ns; ++i) { - p_dKeqdT[i] = 1. - p_dKeqdT[i]; - p_dKeqdT[i] /= (p_group->getT()); - } + 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); } } diff --git a/src/kinetics/RateManager.cpp b/src/kinetics/RateManager.cpp index 557a541e..f11e346e 100644 --- a/src/kinetics/RateManager.cpp +++ b/src/kinetics/RateManager.cpp @@ -127,11 +127,11 @@ RateManager::RateManager(size_t ns, const std::vector& reactions) // 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 + block_size; + mp_dkfdT = mp_gibbs + ns; mp_dkbdT = mp_dkfdT + m_nr; mp_dKeqdT = mp_dkbdT + m_nr; 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/.StateModel.h.swp b/src/thermo/.StateModel.h.swp new file mode 100644 index 0000000000000000000000000000000000000000..af65d2ee722670e5a40b38d023020e594982ea98 GIT binary patch literal 16384 zcmeHNTaO$^6)p!t3`q##B5;J@u#V7PXVx|#yz5vVuh;gGtf!r4>3}J0ftrLx?}7r+a_<_?1K7xG6;9bbC(U;gZ)pnt44v zv*Ng?mAX2!o~T-&TqoGu>}Jtn{fyLCQx(@2&No(Ue$2JRQfjvkriW{1ZUFsSOzQuZ(Ih__-^qTn15F>g{(Lq7&*UM9R1>Yc`obyrzj7L@|Q>2|GOwJ@8h?M zhGqL#it@J??Y=!yUR<^(%YbFTGGH073|Iy%1C{~HfMvikU>UFsSO#7P1CA%e`_b+o zKLEi0|K*MS#+=Yg*QJwO6uz^`$`{tMtq;1j^Tz|Zd&;w9iIpb0z( z{PI2_eg-@Vct8Vq1ULbF5V#9C00`jE_X_b-;3vR0foFg(0Z#x4a1987bHKyE2Y`ow zL%@50@4romU0?$E)mw%5JkSL^;1KY`gF@T@jsqV99suqG4g&wWM~Hs{e+6C!z6*Q{ z_%h%Cp9C%eM}aBe)wc-oW8f>m)4)aGy}*OOZ|@f3N5Hc{7gz?~4gBfNLi`T+4e%21 z9pFXatH3kB<3Jl&03HSo0|$Wby$NFgp94M%JOO+Zm;wF)|NR~K3-BE9MW78hz!t#% zy$876L3~V7qpj1#aZ?Iq+EX5-inL4#?_-eWmoZ3n%#YB2EX@*0J3Zg+k*~>IU)7LF6T8w1U-My_%(iH#S zR^`=Shh2oj7+Fl8f8)U0FA~Iu+i*;zJ2Kf0`HdKT*}qcL5)^v#!c@dZMP$_8*ty6Z14RvytomQ(~4Zgp) z*;BVEzEoYW2KJoPp;*CwnP_9!b@SQxlQeSzx>Y0Po>n5mBEr-^A? zDX1DR&iOUY#&MItxo?J?)o7mE;M~PKUDHL<8!5;c-*bt3+2`9+v6DC9{SJ-h=17a&M2wJs@pdC8^j%WzPf0f|>XuCwSiACogRc3h?*vFI zoy%kuf~Go4C;HK;H#dhye$=IOH*whtdyHN(xX&ASW^fvl#M&chcQF& zPfvCVc3Xz)!hZc>Vh79%5hK!lcTQXvZot}EjVFxQB4fZ_syqc|ud;$;X^94&*e@cNEf@=Hy|JO93*Zx>tNALLm;}WL=NFqRYug*GMr(0DL(z#|n#ObW`p(Wyts7-Em2~US85wD*>nBY_FXunu6@}z8t}M-- znnFk{btqElG&Y+i^R`q*x@R&MZy-x(f-C(#g1~{t;@xr=dTw=-mJwsfq;o6~h~IT1 zaj=ZNcH46J)0k%jDpUwF1SmY4r{BtAliy5kHBK)rFSWK<7S1lURu?xm=#3NBgI@Y6dPT4!YpstCWfLof(I9Mt(!B{v4M_DM#S2>&5V?up?WJ^d_Q*$Oj?NyPr9SSKE;vx$h?cOQ z@KgBMSb8?de2>oLI}lCd7cVSLKtCW0g)Bj$k(X!8;>sQ=DK@%fl0ILWc>+9#Q^x+P z6mQhC%<=vI4{*0@Q|BtWccoDy!1D*!H06YcU06q<5KngU0Ip8?(5Ws!B53q-2 zz%pPNunbrRECZGS%YbFTGGH073|I! 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 +95,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..4f94d3a3 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) @@ -280,7 +284,44 @@ 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 +443,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..77e235d1 100644 --- a/src/transfer/OmegaCE.cpp +++ b/src/transfer/OmegaCE.cpp @@ -44,11 +44,17 @@ class OmegaCE : public TransferModel : 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; + delete [] mp_wrk2; + delete [] mp_wrk3; + delete [] mp_wrk4; } double source() @@ -57,9 +63,34 @@ class OmegaCE : public TransferModel 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(); //Since only omega[0] is taken as source term, jacRho only considers elements of omega[0] - RSCD + } + 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_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..2a6c33a7 100644 --- a/src/transfer/OmegaCElec.cpp +++ b/src/transfer/OmegaCElec.cpp @@ -46,14 +46,24 @@ class OmegaCElec : public TransferModel 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() @@ -69,9 +79,53 @@ class OmegaCElec : public TransferModel 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: + 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..a1a5d62a 100644 --- a/src/transfer/OmegaCV.cpp +++ b/src/transfer/OmegaCV.cpp @@ -60,14 +60,23 @@ class OmegaCV : public TransferModel : 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; + 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$ @@ -88,10 +97,31 @@ class OmegaCV : public TransferModel } } +/** + * 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); + 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(); }; @@ -120,12 +150,73 @@ double const OmegaCV::compute_source_Candler() 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); + for(int i = 0 ; i < m_ns; ++i) + sum += mp_wrk1[i]*mp_wrk2[i]/mixture().speciesMw(i); return(c1*sum*mixture().T()*RU); } +/** + * Assuming non-preferential model for now, density-based jacobians are computed. + * + */ + +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); +} + +/** + * Assuming non-preferential model for now, temperature-based jacobians are computed. + * + */ + +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< OmegaCV, TransferModel> omegaCV("OmegaCV"); diff --git a/src/transfer/OmegaET.cpp b/src/transfer/OmegaET.cpp index 71baf7ec..0e9cc92f 100644 --- a/src/transfer/OmegaET.cpp +++ b/src/transfer/OmegaET.cpp @@ -67,7 +67,7 @@ class OmegaET : public TransferModel if (!m_has_electrons) return 0.0; - const double * p_X = mixture().X(); + const double * p_X = mixture().X(); double nd = mixture().numberDensity(); double T = mixture().T(); double Te = mixture().Te(); @@ -86,6 +86,66 @@ class OmegaET : public TransferModel 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; 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: From c3c0c51cab2fc42f6c24e7d2d37135a131fd5a16 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Fri, 6 Jun 2025 02:17:16 -0700 Subject: [PATCH 14/23] Remove unnecessary files from tracking and update .gitignore --- .gitignore | 4 ++++ compile.sh | 4 ---- data/mechanisms/air5p_Park.xml | 29 ----------------------------- 3 files changed, 4 insertions(+), 33 deletions(-) delete mode 100755 compile.sh delete mode 100644 data/mechanisms/air5p_Park.xml diff --git a/.gitignore b/.gitignore index 7e372b97..7e59278b 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,7 @@ _skbuild/ # OS specifics .DS_Store + +#Files to ignore +echo "compile.sh" >> .gitignore +echo "data/mechanisms/air5p_Park.xml" >> .gitignore diff --git a/compile.sh b/compile.sh deleted file mode 100755 index 9a1bbad6..00000000 --- a/compile.sh +++ /dev/null @@ -1,4 +0,0 @@ -#!/bin/bash -rm -r build; mkdir build; cd build -cmake -DCMAKE_INSTALL_PREFIX=/Users/raghava/Desktop/CHANL_HPC/Mutationpp/install .. -make -j 20 install diff --git a/data/mechanisms/air5p_Park.xml b/data/mechanisms/air5p_Park.xml deleted file mode 100644 index 584e7a84..00000000 --- a/data/mechanisms/air5p_Park.xml +++ /dev/null @@ -1,29 +0,0 @@ - - - - - - - - - N2:0.2333, NO:0.2333, O2:0.2333, N:1.0, O:1.0 - - - - - - - - - - - - - - - From 9ec353f9128f12159dbb2fad2564efa2680a1b97 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Wed, 30 Jul 2025 02:32:53 -0700 Subject: [PATCH 15/23] Added argon.xml to data/mixtures --- data/mixtures/argon.xml | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 data/mixtures/argon.xml diff --git a/data/mixtures/argon.xml b/data/mixtures/argon.xml new file mode 100644 index 00000000..846bbe51 --- /dev/null +++ b/data/mixtures/argon.xml @@ -0,0 +1,11 @@ + + + + Ar + + + + Ar: 1.0 + + + From cfedb481da5f598304a833922e220dfe2b3936c6 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Wed, 13 Aug 2025 04:33:45 -0700 Subject: [PATCH 16/23] Modified CollisionDB.cpp to avoid Nan errors when mixture is a single species. Allows to run SU2 to run for single species. Also, added nitrogen (N2, N) mixture and N2 dissociation mechanism in data. --- data/mechanisms/nitrogen_dissociation.xml | 11 +++++++++++ data/mixtures/nitrogen.xml | 13 +++++++++++++ src/transport/CollisionDB.cpp | 4 ++++ 3 files changed, 28 insertions(+) create mode 100644 data/mechanisms/nitrogen_dissociation.xml create mode 100644 data/mixtures/nitrogen.xml diff --git a/data/mechanisms/nitrogen_dissociation.xml b/data/mechanisms/nitrogen_dissociation.xml new file mode 100644 index 00000000..cdae207b --- /dev/null +++ b/data/mechanisms/nitrogen_dissociation.xml @@ -0,0 +1,11 @@ + + + + + + + + N2:0.2333 + + + diff --git a/data/mixtures/nitrogen.xml b/data/mixtures/nitrogen.xml new file mode 100644 index 00000000..dfe2bf1c --- /dev/null +++ b/data/mixtures/nitrogen.xml @@ -0,0 +1,13 @@ + + + + + N N2 + + + + N:1.0 + + + + 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) { From cb11b883d583e8a0ef74ce74376d3b29a4bb1432 Mon Sep 17 00:00:00 2001 From: raghava-davuluri Date: Mon, 18 Aug 2025 19:54:25 -0700 Subject: [PATCH 17/23] Update .gitignore --- .gitignore | 4 ---- 1 file changed, 4 deletions(-) diff --git a/.gitignore b/.gitignore index 7e59278b..7e372b97 100644 --- a/.gitignore +++ b/.gitignore @@ -21,7 +21,3 @@ _skbuild/ # OS specifics .DS_Store - -#Files to ignore -echo "compile.sh" >> .gitignore -echo "data/mechanisms/air5p_Park.xml" >> .gitignore From 50678e79ec41c8acff172482a456c3ff96af7690 Mon Sep 17 00:00:00 2001 From: raghava-davuluri Date: Mon, 18 Aug 2025 20:01:20 -0700 Subject: [PATCH 18/23] Delete src/thermo/.StateModel.h.swp --- src/thermo/.StateModel.h.swp | Bin 16384 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/thermo/.StateModel.h.swp diff --git a/src/thermo/.StateModel.h.swp b/src/thermo/.StateModel.h.swp deleted file mode 100644 index af65d2ee722670e5a40b38d023020e594982ea98..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16384 zcmeHNTaO$^6)p!t3`q##B5;J@u#V7PXVx|#yz5vVuh;gGtf!r4>3}J0ftrLx?}7r+a_<_?1K7xG6;9bbC(U;gZ)pnt44v zv*Ng?mAX2!o~T-&TqoGu>}Jtn{fyLCQx(@2&No(Ue$2JRQfjvkriW{1ZUFsSOzQuZ(Ih__-^qTn15F>g{(Lq7&*UM9R1>Yc`obyrzj7L@|Q>2|GOwJ@8h?M zhGqL#it@J??Y=!yUR<^(%YbFTGGH073|Iy%1C{~HfMvikU>UFsSO#7P1CA%e`_b+o zKLEi0|K*MS#+=Yg*QJwO6uz^`$`{tMtq;1j^Tz|Zd&;w9iIpb0z( z{PI2_eg-@Vct8Vq1ULbF5V#9C00`jE_X_b-;3vR0foFg(0Z#x4a1987bHKyE2Y`ow zL%@50@4romU0?$E)mw%5JkSL^;1KY`gF@T@jsqV99suqG4g&wWM~Hs{e+6C!z6*Q{ z_%h%Cp9C%eM}aBe)wc-oW8f>m)4)aGy}*OOZ|@f3N5Hc{7gz?~4gBfNLi`T+4e%21 z9pFXatH3kB<3Jl&03HSo0|$Wby$NFgp94M%JOO+Zm;wF)|NR~K3-BE9MW78hz!t#% zy$876L3~V7qpj1#aZ?Iq+EX5-inL4#?_-eWmoZ3n%#YB2EX@*0J3Zg+k*~>IU)7LF6T8w1U-My_%(iH#S zR^`=Shh2oj7+Fl8f8)U0FA~Iu+i*;zJ2Kf0`HdKT*}qcL5)^v#!c@dZMP$_8*ty6Z14RvytomQ(~4Zgp) z*;BVEzEoYW2KJoPp;*CwnP_9!b@SQxlQeSzx>Y0Po>n5mBEr-^A? zDX1DR&iOUY#&MItxo?J?)o7mE;M~PKUDHL<8!5;c-*bt3+2`9+v6DC9{SJ-h=17a&M2wJs@pdC8^j%WzPf0f|>XuCwSiACogRc3h?*vFI zoy%kuf~Go4C;HK;H#dhye$=IOH*whtdyHN(xX&ASW^fvl#M&chcQF& zPfvCVc3Xz)!hZc>Vh79%5hK!lcTQXvZot}EjVFxQB4fZ_syqc|ud;$;X^94&*e@cNEf@=Hy|JO93*Zx>tNALLm;}WL=NFqRYug*GMr(0DL(z#|n#ObW`p(Wyts7-Em2~US85wD*>nBY_FXunu6@}z8t}M-- znnFk{btqElG&Y+i^R`q*x@R&MZy-x(f-C(#g1~{t;@xr=dTw=-mJwsfq;o6~h~IT1 zaj=ZNcH46J)0k%jDpUwF1SmY4r{BtAliy5kHBK)rFSWK<7S1lURu?xm=#3NBgI@Y6dPT4!YpstCWfLof(I9Mt(!B{v4M_DM#S2>&5V?up?WJ^d_Q*$Oj?NyPr9SSKE;vx$h?cOQ z@KgBMSb8?de2>oLI}lCd7cVSLKtCW0g)Bj$k(X!8;>sQ=DK@%fl0ILWc>+9#Q^x+P z6mQhC%<=vI4{*0@Q|BtWccoDy!1D*!H06YcU06q<5KngU0Ip8?(5Ws!B53q-2 zz%pPNunbrRECZGS%YbFTGGH073|I! Date: Mon, 18 Aug 2025 20:07:42 -0700 Subject: [PATCH 19/23] Resolved issues raised by reviewers --- .gitignore | 4 ---- data/mechanisms/nitrogen_dissociation.xml | 11 ----------- data/mixtures/argon.xml | 11 ----------- data/mixtures/nitrogen.xml | 13 ------------- src/thermo/.StateModel.h.swp | Bin 16384 -> 0 bytes src/transfer/OmegaCE.cpp | 6 +++--- 6 files changed, 3 insertions(+), 42 deletions(-) delete mode 100644 data/mechanisms/nitrogen_dissociation.xml delete mode 100644 data/mixtures/argon.xml delete mode 100644 data/mixtures/nitrogen.xml delete mode 100644 src/thermo/.StateModel.h.swp diff --git a/.gitignore b/.gitignore index 7e59278b..7e372b97 100644 --- a/.gitignore +++ b/.gitignore @@ -21,7 +21,3 @@ _skbuild/ # OS specifics .DS_Store - -#Files to ignore -echo "compile.sh" >> .gitignore -echo "data/mechanisms/air5p_Park.xml" >> .gitignore diff --git a/data/mechanisms/nitrogen_dissociation.xml b/data/mechanisms/nitrogen_dissociation.xml deleted file mode 100644 index cdae207b..00000000 --- a/data/mechanisms/nitrogen_dissociation.xml +++ /dev/null @@ -1,11 +0,0 @@ - - - - - - - - N2:0.2333 - - - diff --git a/data/mixtures/argon.xml b/data/mixtures/argon.xml deleted file mode 100644 index 846bbe51..00000000 --- a/data/mixtures/argon.xml +++ /dev/null @@ -1,11 +0,0 @@ - - - - Ar - - - - Ar: 1.0 - - - diff --git a/data/mixtures/nitrogen.xml b/data/mixtures/nitrogen.xml deleted file mode 100644 index dfe2bf1c..00000000 --- a/data/mixtures/nitrogen.xml +++ /dev/null @@ -1,13 +0,0 @@ - - - - - N N2 - - - - N:1.0 - - - - diff --git a/src/thermo/.StateModel.h.swp b/src/thermo/.StateModel.h.swp deleted file mode 100644 index af65d2ee722670e5a40b38d023020e594982ea98..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16384 zcmeHNTaO$^6)p!t3`q##B5;J@u#V7PXVx|#yz5vVuh;gGtf!r4>3}J0ftrLx?}7r+a_<_?1K7xG6;9bbC(U;gZ)pnt44v zv*Ng?mAX2!o~T-&TqoGu>}Jtn{fyLCQx(@2&No(Ue$2JRQfjvkriW{1ZUFsSOzQuZ(Ih__-^qTn15F>g{(Lq7&*UM9R1>Yc`obyrzj7L@|Q>2|GOwJ@8h?M zhGqL#it@J??Y=!yUR<^(%YbFTGGH073|Iy%1C{~HfMvikU>UFsSO#7P1CA%e`_b+o zKLEi0|K*MS#+=Yg*QJwO6uz^`$`{tMtq;1j^Tz|Zd&;w9iIpb0z( z{PI2_eg-@Vct8Vq1ULbF5V#9C00`jE_X_b-;3vR0foFg(0Z#x4a1987bHKyE2Y`ow zL%@50@4romU0?$E)mw%5JkSL^;1KY`gF@T@jsqV99suqG4g&wWM~Hs{e+6C!z6*Q{ z_%h%Cp9C%eM}aBe)wc-oW8f>m)4)aGy}*OOZ|@f3N5Hc{7gz?~4gBfNLi`T+4e%21 z9pFXatH3kB<3Jl&03HSo0|$Wby$NFgp94M%JOO+Zm;wF)|NR~K3-BE9MW78hz!t#% zy$876L3~V7qpj1#aZ?Iq+EX5-inL4#?_-eWmoZ3n%#YB2EX@*0J3Zg+k*~>IU)7LF6T8w1U-My_%(iH#S zR^`=Shh2oj7+Fl8f8)U0FA~Iu+i*;zJ2Kf0`HdKT*}qcL5)^v#!c@dZMP$_8*ty6Z14RvytomQ(~4Zgp) z*;BVEzEoYW2KJoPp;*CwnP_9!b@SQxlQeSzx>Y0Po>n5mBEr-^A? zDX1DR&iOUY#&MItxo?J?)o7mE;M~PKUDHL<8!5;c-*bt3+2`9+v6DC9{SJ-h=17a&M2wJs@pdC8^j%WzPf0f|>XuCwSiACogRc3h?*vFI zoy%kuf~Go4C;HK;H#dhye$=IOH*whtdyHN(xX&ASW^fvl#M&chcQF& zPfvCVc3Xz)!hZc>Vh79%5hK!lcTQXvZot}EjVFxQB4fZ_syqc|ud;$;X^94&*e@cNEf@=Hy|JO93*Zx>tNALLm;}WL=NFqRYug*GMr(0DL(z#|n#ObW`p(Wyts7-Em2~US85wD*>nBY_FXunu6@}z8t}M-- znnFk{btqElG&Y+i^R`q*x@R&MZy-x(f-C(#g1~{t;@xr=dTw=-mJwsfq;o6~h~IT1 zaj=ZNcH46J)0k%jDpUwF1SmY4r{BtAliy5kHBK)rFSWK<7S1lURu?xm=#3NBgI@Y6dPT4!YpstCWfLof(I9Mt(!B{v4M_DM#S2>&5V?up?WJ^d_Q*$Oj?NyPr9SSKE;vx$h?cOQ z@KgBMSb8?de2>oLI}lCd7cVSLKtCW0g)Bj$k(X!8;>sQ=DK@%fl0ILWc>+9#Q^x+P z6mQhC%<=vI4{*0@Q|BtWccoDy!1D*!H06YcU06q<5KngU0Ip8?(5Ws!B53q-2 zz%pPNunbrRECZGS%YbFTGGH073|I! Date: Wed, 8 Oct 2025 21:10:42 -0700 Subject: [PATCH 20/23] Modified indentations from tabs to spaces --- src/kinetics/JacobianManager.cpp | 324 +++++++++++++++---------------- src/thermo/StateModel.h | 42 ++-- src/transfer/OmegaCE.cpp | 93 ++++----- src/transfer/OmegaCElec.cpp | 138 ++++++------- src/transfer/OmegaCV.cpp | 184 +++++++++--------- src/transfer/OmegaET.cpp | 84 ++++---- 6 files changed, 428 insertions(+), 437 deletions(-) diff --git a/src/kinetics/JacobianManager.cpp b/src/kinetics/JacobianManager.cpp index aa770925..05932d88 100644 --- a/src/kinetics/JacobianManager.cpp +++ b/src/kinetics/JacobianManager.cpp @@ -347,150 +347,148 @@ 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; - } - } - -} + { + + 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; - } - } - -} + { + + 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; + } + } + } //============================================================================== @@ -546,38 +544,36 @@ void JacobianManager::computeJacobian( 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); - - computeTReactionTDerivatives(reactions, mp_T, dTrxndT); - - for (int i = 0; i < reactions.size(); ++i) - p_dropT[i] = p_dropf[i]*dTrxndT[i].first - p_dropb[i]*dTrxndT[i].second; - -} + { + const size_t nt = m_thermo.nEnergyEqns(); + std::vector > dTrxndT; + + std::fill(mp_T, mp_T + nt, 0.); + m_thermo.getTemperatures(mp_T); + + computeTReactionTDerivatives(reactions, mp_T, dTrxndT); + + 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); - - computeTReactionTvDerivatives(reactions, mp_T, dTrxndTv); - - for (int i = 0; i < reactions.size(); ++i) - p_dropTv[i] = p_dropf[i]*dTrxndTv[i].first - p_dropb[i]*dTrxndTv[i].second; - -} + { + const size_t nt = m_thermo.nEnergyEqns(); + std::vector > dTrxndTv; + + std::fill(mp_T, mp_T + nt, 0.); + m_thermo.getTemperatures(mp_T); + + computeTReactionTvDerivatives(reactions, mp_T, dTrxndTv); + + for (int i = 0; i < reactions.size(); ++i) + p_dropTv[i] = p_dropf[i]*dTrxndTv[i].first - p_dropb[i]*dTrxndTv[i].second; + } //============================================================================== diff --git a/src/thermo/StateModel.h b/src/thermo/StateModel.h index 4f94d3a3..d213cfb2 100644 --- a/src/thermo/StateModel.h +++ b/src/thermo/StateModel.h @@ -214,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()"); @@ -223,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()"); @@ -232,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){ @@ -242,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()"); @@ -292,16 +292,15 @@ class StateModel 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; - + 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]; + p_omegaJRho[j] += mp_omegaRho[j]; } - } /** @@ -310,16 +309,15 @@ class StateModel */ 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; - + 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]; + p_omegaJTTv[j] += mp_omegaTTv[j]; } - } protected: diff --git a/src/transfer/OmegaCE.cpp b/src/transfer/OmegaCE.cpp index c77fe06d..c40aca4a 100644 --- a/src/transfer/OmegaCE.cpp +++ b/src/transfer/OmegaCE.cpp @@ -40,57 +40,58 @@ namespace Mutation { class OmegaCE : public TransferModel { public: - 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(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; - delete [] mp_wrk2; - delete [] mp_wrk3; - delete [] mp_wrk4; - } + ~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(); //Since only omega[0] is taken as source term, jacRho only considers elements of omega[0] - RSCD - } - 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]; - } + 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_wrk2; - double* mp_wrk3; - double* mp_wrk4; + 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 2a6c33a7..1d20d9c1 100644 --- a/src/transfer/OmegaCElec.cpp +++ b/src/transfer/OmegaCElec.cpp @@ -43,89 +43,89 @@ namespace Mutation { class OmegaCElec : public TransferModel { public: - 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); + 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.); + 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); + mixture().speciesHOverRT(NULL, NULL, NULL, NULL, mp_wrk1, NULL); + mixture().jacobianRho(mp_wrk3); - double aux = mixture().T()*RU; + 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); + 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.); + void jacobianTTv(double* const p_jacTTv) + { + std::fill(p_jacTTv, p_jacTTv + m_nt, 0.); - double T=mixture().T(); - double Tv=mixture().Tv(); + 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); + 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; + 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); - } - + 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: - int m_ns; - int m_nt; - double* mp_wrk1; - double* mp_wrk2; - double* mp_wrk3; - double* mp_wrk4; - double* mp_wrk5; - double* mp_wrk6; + 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 a1a5d62a..f18fba12 100644 --- a/src/transfer/OmegaCV.cpp +++ b/src/transfer/OmegaCV.cpp @@ -56,46 +56,46 @@ class OmegaCV : public TransferModel */ public: - 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(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; - delete [] mp_wrk3; - delete [] mp_wrk4; - delete [] mp_wrk5; - delete [] mp_wrk6; - }; + ~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; + } + } /** * Computes the density-based jacobian terms of the Vibration-Chemistry energy transfer.$ @@ -103,7 +103,7 @@ class OmegaCV : public TransferModel * \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); + void jacobianRho(double* const p_jacRho); /** * Computes the temperature-based jacobian terms of the Vibration-Chemistry energy transfer.$ @@ -111,19 +111,19 @@ class OmegaCV : public TransferModel * \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); + void jacobianTTv(double* const p_jacTTv); 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(); + 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(); }; /** @@ -136,86 +136,86 @@ class OmegaCV : public TransferModel * */ -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); - } + 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); + } /** * Assuming non-preferential model for now, density-based jacobians are computed. * */ -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); - + 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); - + 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 j = 0; j < m_ns; ++j) + p_jacRho[i] += aux*mp_wrk1[i]*mp_wrk3[j*m_ns+i]/mixture().speciesMw(i); + } /** * Assuming non-preferential model for now, temperature-based jacobians are computed. * */ -void OmegaCV::jacobianTTv(double* const p_jacTTv) -{ - std::fill(p_jacTTv, p_jacTTv + m_nt, 0.); - + 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 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); + 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 0e9cc92f..bcdebe68 100644 --- a/src/transfer/OmegaET.cpp +++ b/src/transfer/OmegaET.cpp @@ -57,98 +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) - { + } + 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(); + 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; - + + 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) - { - + 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(); + 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; }; From 87002b41b3cc682aa08e7aa02a5ae3f3f152be60 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Thu, 9 Oct 2025 13:03:20 -0700 Subject: [PATCH 21/23] Remove some brackets --- src/kinetics/RateManager.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/kinetics/RateManager.cpp b/src/kinetics/RateManager.cpp index f11e346e..c778aa3c 100644 --- a/src/kinetics/RateManager.cpp +++ b/src/kinetics/RateManager.cpp @@ -121,9 +121,8 @@ RateManager::RateManager(size_t ns, const std::vector& reactions) { // 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) { + 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 From 2ac8195ff12496d7133a648b4fe7dbe59d9c4dc6 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Sat, 11 Oct 2025 00:24:47 -0700 Subject: [PATCH 22/23] Modified tests to include chemistry jacobian tests --- tests/c++/test_reaction_mechanism.cpp | 126 +++++++++++++++++++++++++- tests/c++/test_reaction_rates.cpp | 12 +-- 2 files changed, 131 insertions(+), 7 deletions(-) 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 From 58a82084c710b8b9ea58d4276f3a39d2270bb61d Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Fri, 12 Dec 2025 15:36:46 -0500 Subject: [PATCH 23/23] Modified Kinetics to give proper jacobians for 1T model and modified some lines based on space restrictions. --- src/kinetics/JacobianManager.cpp | 18 ++++++++++++++---- src/kinetics/Kinetics.cpp | 3 ++- src/thermo/HarmonicOscillator.cpp | 4 +++- src/transfer/OmegaCV.cpp | 8 +++++--- 4 files changed, 24 insertions(+), 9 deletions(-) diff --git a/src/kinetics/JacobianManager.cpp b/src/kinetics/JacobianManager.cpp index 05932d88..bb716414 100644 --- a/src/kinetics/JacobianManager.cpp +++ b/src/kinetics/JacobianManager.cpp @@ -550,9 +550,14 @@ void JacobianManager::computeJacobianT( std::fill(mp_T, mp_T + nt, 0.); m_thermo.getTemperatures(mp_T); - - computeTReactionTDerivatives(reactions, mp_T, dTrxndT); - + + 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; } @@ -569,7 +574,12 @@ void JacobianManager::computeJacobianTv( std::fill(mp_T, mp_T + nt, 0.); m_thermo.getTemperatures(mp_T); - computeTReactionTvDerivatives(reactions, mp_T, dTrxndTv); + 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; diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 04b46926..b5d17b63 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -306,7 +306,7 @@ void Kinetics::forwardRateOfProgressDerivatives(double* const p_dropf) return; mp_rates->update(m_thermo); - ArrayXd dkfdT = + ArrayXd dkfdT = Map(mp_rates->dkfdT(), nReactions()); forwardRatesOfProgress(mp_ropf); @@ -498,6 +498,7 @@ void Kinetics::jacobianRho(double* const p_jac) // Compute the Jacobian matrix m_jacobian.computeJacobian(mp_ropf, mp_ropb, mp_rop, p_jac); + } //============================================================================== diff --git a/src/thermo/HarmonicOscillator.cpp b/src/thermo/HarmonicOscillator.cpp index 7bcba079..de8294d7 100644 --- a/src/thermo/HarmonicOscillator.cpp +++ b/src/thermo/HarmonicOscillator.cpp @@ -44,7 +44,9 @@ 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); + energyTderivative += + std::pow(theta / T, 2) * std::exp(theta / T) / + std::pow(std::exp(theta / T) - 1.0, 2); return energyTderivative; } diff --git a/src/transfer/OmegaCV.cpp b/src/transfer/OmegaCV.cpp index f18fba12..bd1089d8 100644 --- a/src/transfer/OmegaCV.cpp +++ b/src/transfer/OmegaCV.cpp @@ -100,7 +100,8 @@ class OmegaCV : public TransferModel /** * 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] + * \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); @@ -108,8 +109,9 @@ class OmegaCV : public TransferModel /** * 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] + * \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);