Skip to content

Analytical jacobians for chemistry and vibrational source terms #289

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
raghava-davuluri wants to merge 24 commits into
base: master
Choose a base branch
from
Open

Analytical jacobians for chemistry and vibrational source terms #289

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
87e1c1b
New file from Catarina that helps with NAN in SU2 coupling, See 1/20/…
Jan 20, 2025
e436e0c
Successfully builds; added function arguments to setState()
aaronmlarsen Jan 21, 2025
382b64d
Updated function headers to include Tv_old and T_old to protect again…
aaronmlarsen Jan 21, 2025
bd2f482
Fixed redefinition of default arguments
aaronmlarsen Jan 21, 2025
2b41465
Builds successfully. Removes T_old, Tv_old
aaronmlarsen Jan 21, 2025
e8e719d
Revert "Successfully builds; added function arguments to setState()"
aaronmlarsen Feb 3, 2025
58e4e40
Revert "Builds successfully. Removes T_old, Tv_old"
aaronmlarsen Feb 3, 2025
c08f303
Revert "Fixed redefinition of default arguments"
aaronmlarsen Feb 3, 2025
b43bbb0
Revert "Updated function headers to include Tv_old and T_old to prote…
aaronmlarsen Feb 3, 2025
e3b5c19
Revert "New file from Catarina that helps with NAN in SU2 coupling, S…
aaronmlarsen Feb 3, 2025
405b6d4
Added some variables and functions to caluclate the temperature-depen…
Apr 9, 2025
c307612
Finished adding necessary functions and variables that will help in c…
Apr 16, 2025
dac2a5d
Added functions and variables that help in computing the vibrational …
raghava-davuluri Jun 6, 2025
c3c0c51
Remove unnecessary files from tracking and update .gitignore
raghava-davuluri Jun 6, 2025
9ec353f
Added argon.xml to data/mixtures
raghava-davuluri Jul 30, 2025
cfedb48
Modified CollisionDB.cpp to avoid Nan errors when mixture is a single…
raghava-davuluri Aug 13, 2025
cb11b88
Update .gitignore
raghava-davuluri Aug 19, 2025
50678e7
Delete src/thermo/.StateModel.h.swp
raghava-davuluri Aug 19, 2025
bd89749
Resolved issues raised by reviewers
raghava-davuluri Aug 19, 2025
aa0bb5a
Merge branch 'feature_jacobians' of github.com:hypersonic-lab/Mutatio…
raghava-davuluri Aug 19, 2025
9d07fa1
Modified indentations from tabs to spaces
raghava-davuluri Oct 9, 2025
87002b4
Remove some brackets
raghava-davuluri Oct 9, 2025
2ac8195
Modified tests to include chemistry jacobian tests
raghava-davuluri Oct 11, 2025
58a8208
Modified Kinetics to give proper jacobians for 1T model and modified …
raghava-davuluri Dec 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion data/transport/collisions.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 </Q22>
</pair>

</collisions>
</collisions>
16 changes: 16 additions & 0 deletions src/general/Mixture.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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().
Expand Down
276 changes: 273 additions & 3 deletions src/kinetics/JacobianManager.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,29 @@ void ReactionStoich<Reactants, Products>::contributeToJacobian(

//==============================================================================

template <typename Reactants, typename Products>
void ReactionStoich<Reactants, Products>::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 <typename Reactants, typename Products>
void ThirdbodyReactionStoich<Reactants, Products>::contributeToJacobian(
const double kf, const double kb, const double* const conc,
Expand All @@ -104,6 +127,34 @@ void ThirdbodyReactionStoich<Reactants, Products>::contributeToJacobian(

//==============================================================================

template <typename Reactants, typename Products>
void ThirdbodyReactionStoich<Reactants, Products>::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 <typename Reactants>
void JacobianManager::addReactionStoich(
JacStoichBase* p_reacs, JacStoichBase* p_prods, const StoichType type,
Expand Down Expand Up @@ -293,6 +344,179 @@ void JacobianManager::addReaction(const Reaction& reaction)

//==============================================================================

void JacobianManager::computeTReactionTDerivatives(
const std::vector<Reaction>& reactions, const double* const p_t,
std::vector<std::pair<double, double> >& 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<Reaction>& reactions, const double* const p_t,
std::vector<std::pair<double, double> >& dTrxndTv) const
{

for(size_t i = 0; i < reactions.size(); ++i) {
const size_t type = reactions[i].type();
switch (type) {
case ASSOCIATIVE_IONIZATION:
dTrxndTv.push_back({0.0, 1.0});
break;
case DISSOCIATIVE_RECOMBINATION:
dTrxndTv.push_back({1.0, 0.0});
break;
case ASSOCIATIVE_DETACHMENT:
dTrxndTv.push_back({0.0, 1.0});
break;
case DISSOCIATIVE_ATTACHMENT:
dTrxndTv.push_back({1.0, 0.0});
break;
case DISSOCIATION_M:
dTrxndTv.push_back({0.5*sqrt(p_t[0]*p_t[1])/p_t[1], 0.0});
break;
case DISSOCIATION_E:
dTrxndTv.push_back({1.0, 1.0});
break;
case RECOMBINATION_M:
dTrxndTv.push_back({0.0, 0.5*sqrt(p_t[0]*p_t[1])/p_t[1]});
break;
case RECOMBINATION_E:
dTrxndTv.push_back({1.0, 1.0});
break;
case IONIZATION_M:
dTrxndTv.push_back({0.0, 0.0});
break;
case IONIZATION_E:
dTrxndTv.push_back({1.0, 0.0});
break;
case ION_RECOMBINATION_M:
dTrxndTv.push_back({0.0, 0.0});
break;
case ION_RECOMBINATION_E:
dTrxndTv.push_back({0.0, 1.0});
break;
case ELECTRONIC_ATTACHMENT_E:
dTrxndTv.push_back({1.0, 1.0});
break;
case ELECTRONIC_DETACHMENT_E:
dTrxndTv.push_back({1.0, 1.0});
break;
case ELECTRONIC_ATTACHMENT_M:
dTrxndTv.push_back({1.0, 0.0});
break;
case ELECTRONIC_DETACHMENT_M:
dTrxndTv.push_back({0.0, 1.0});
break;
case EXCHANGE:
dTrxndTv.push_back({0.0, 0.0});
break;
case EXCITATION_E:
dTrxndTv.push_back({1.0, 1.0});
break;
case EXCITATION_M:
dTrxndTv.push_back({0.0, 0.0});
break;
default:
dTrxndTv.push_back({0.0, 0.0});
break;
}
}
}

//==============================================================================

void JacobianManager::computeNetRateProgressJacobian(
const double* const kf, const double* const kb, const double* const conc,
double* const p_dropRho) const
{
const size_t ns = m_thermo.nSpecies();
const size_t nr = m_reactions.size();

// Make sure we are staring with a clean slate
std::fill(p_dropRho, p_dropRho + ns*nr, 0.0);

// Loop over each reaction and compute the dRR/dconc_k
for (int i = 0; i < nr; ++i)
m_reactions[i]->contributeToNetProgressJacobian(
kf[i], kb[i], conc, mp_work1, p_dropRho, ns, i);


// Finally, multiply by the species molecular weight ratios
for (int i = 0, index = 0; i < nr; ++i)
for (int j = 0; j < ns; ++j, ++index)
p_dropRho[index] /= m_thermo.speciesMw(j);

}

//==============================================================================

void JacobianManager::computeJacobian(
const double* const kf, const double* const kb, const double* const conc,
double* const sjac) const
Expand All @@ -304,17 +528,63 @@ void JacobianManager::computeJacobian(
std::fill(sjac, sjac+ns*ns, 0.0);

// Loop over each reaction and compute the dRR/dconc_k
for (int i = 0; i < nr; ++i)
for (int i = 0; i < nr; ++i)
m_reactions[i]->contributeToJacobian(
kf[i], kb[i], conc, mp_work, sjac, ns);


// Finally, multiply by the species molecular weight ratios
for (int i = 0, index = 0; i < ns; ++i) {
for (int i = 0, index = 0; i < ns; ++i)
for (int j = 0; j < ns; ++j, ++index)
sjac[index] *= m_thermo.speciesMw(i) / m_thermo.speciesMw(j);
}
}

//==============================================================================

void JacobianManager::computeJacobianT(
const double* const p_dropf, const double* const p_dropb,
const std::vector<Reaction>& reactions, double* const p_dropT) const
{
const size_t nt = m_thermo.nEnergyEqns();
std::vector<std::pair<double, double> > dTrxndT;

std::fill(mp_T, mp_T + nt, 0.);
m_thermo.getTemperatures(mp_T);

if (nt == 2) {
computeTReactionTDerivatives(reactions, mp_T, dTrxndT);
} else {
for (int i = 0; i < reactions.size(); ++i)
dTrxndT.push_back({1.0, 1.0});
}

for (int i = 0; i < reactions.size(); ++i)
p_dropT[i] = p_dropf[i]*dTrxndT[i].first - p_dropb[i]*dTrxndT[i].second;
}

//==============================================================================

void JacobianManager::computeJacobianTv(
const double* const p_dropf, const double* const p_dropb,
const std::vector<Reaction>& reactions, double* const p_dropTv) const
{
const size_t nt = m_thermo.nEnergyEqns();
std::vector<std::pair<double, double> > dTrxndTv;

std::fill(mp_T, mp_T + nt, 0.);
m_thermo.getTemperatures(mp_T);

if (nt == 2) {
computeTReactionTvDerivatives(reactions, mp_T, dTrxndTv);
} else {
for (int i = 0; i < reactions.size(); ++i)
dTrxndTv.push_back({0.0, 0.0});
}

for (int i = 0; i < reactions.size(); ++i)
p_dropTv[i] = p_dropf[i]*dTrxndTv[i].first - p_dropb[i]*dTrxndTv[i].second;
}

//==============================================================================

} // namespace Kinetics
Expand Down
Loading