/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see .
\*---------------------------------------------------------------------------*/
#include "TDACChemistryModel.H"
#include "UniformField.H"
#include "localEulerDdtScheme.H"
#include "clockTime.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template
Foam::TDACChemistryModel::TDACChemistryModel
(
const ReactionThermo& thermo
)
:
StandardChemistryModel(thermo),
variableTimeStep_
(
this->mesh().time().controlDict().lookupOrDefault
(
"adjustTimeStep",
false
)
|| fv::localEulerDdt::enabled(this->mesh())
),
timeSteps_(0),
NsDAC_(this->nSpecie_),
completeC_(this->nSpecie_, 0),
reactionsDisabled_(this->reactions_.size(), false),
specieComp_(this->nSpecie_),
completeToSimplifiedIndex_(this->nSpecie_, -1),
simplifiedToCompleteIndex_(this->nSpecie_),
tabulationResults_
(
IOobject
(
thermo.phasePropertyName("TabulationResults"),
this->time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
this->mesh(),
scalar(0)
)
{
const basicSpecieMixture& composition = this->thermo().composition();
const HashTable>& specComp =
dynamicCast&>(this->thermo())
.specieComposition();
forAll(specieComp_, i)
{
specieComp_[i] = specComp[this->Y()[i].member()];
}
mechRed_ = chemistryReductionMethod::New
(
*this,
*this
);
// When the mechanism reduction method is used, the 'active' flag for every
// species should be initialized (by default 'active' is true)
if (mechRed_->active())
{
forAll(this->Y(), i)
{
IOobject header
(
this->Y()[i].name(),
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ
);
// Check if the species file is provided, if not set inactive
// and NO_WRITE
if (!header.typeHeaderOk(true))
{
composition.setInactive(i);
}
}
}
tabulation_ = chemistryTabulationMethod::New
(
*this,
*this
);
if (mechRed_->log())
{
cpuReduceFile_ = logFile("cpu_reduce.out");
nActiveSpeciesFile_ = logFile("nActiveSpecies.out");
}
if (tabulation_->log())
{
cpuAddFile_ = logFile("cpu_add.out");
cpuGrowFile_ = logFile("cpu_grow.out");
cpuRetrieveFile_ = logFile("cpu_retrieve.out");
}
if (mechRed_->log() || tabulation_->log())
{
cpuSolveFile_ = logFile("cpu_solve.out");
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template
Foam::TDACChemistryModel::~TDACChemistryModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template
void Foam::TDACChemistryModel::omega
(
const scalar p,
const scalar T,
const scalarField& c, // Contains all species even when mechRed is active
const label li,
scalarField& dcdt
) const
{
const bool reduced = mechRed_->active();
scalar pf, cf, pr, cr;
label lRef, rRef;
dcdt = Zero;
forAll(this->reactions_, i)
{
if (!reactionsDisabled_[i])
{
const Reaction& R = this->reactions_[i];
scalar omegai = R.omega
(
p, T, c, li, pf, cf, lRef, pr, cr, rRef
);
forAll(R.lhs(), s)
{
label si = R.lhs()[s].index;
if (reduced)
{
si = completeToSimplifiedIndex_[si];
}
const scalar sl = R.lhs()[s].stoichCoeff;
dcdt[si] -= sl*omegai;
}
forAll(R.rhs(), s)
{
label si = R.rhs()[s].index;
if (reduced)
{
si = completeToSimplifiedIndex_[si];
}
const scalar sr = R.rhs()[s].stoichCoeff;
dcdt[si] += sr*omegai;
}
}
}
}
template
Foam::scalar Foam::TDACChemistryModel::omega
(
const Reaction& R,
const scalar p,
const scalar T,
const scalarField& c, // Contains all species even when mechRed is active
const label li,
scalar& pf,
scalar& cf,
label& lRef,
scalar& pr,
scalar& cr,
label& rRef
) const
{
const scalar kf = R.kf(p, T, c, li);
const scalar kr = R.kr(kf, p, T, c, li);
const label Nl = R.lhs().size();
const label Nr = R.rhs().size();
label slRef = 0;
lRef = R.lhs()[slRef].index;
pf = kf;
for (label s=1; s small)
{
pf *= pow(cf, exp - 1);
}
else
{
pf = 0;
}
}
else
{
pf *= pow(cf, exp - 1);
}
}
label srRef = 0;
rRef = R.rhs()[srRef].index;
// Find the matrix element and element position for the rhs
pr = kr;
for (label s=1; s small)
{
pr *= pow(cr, exp - 1);
}
else
{
pr = 0;
}
}
else
{
pr *= pow(cr, exp - 1);
}
}
return pf*cf - pr*cr;
}
template
void Foam::TDACChemistryModel::derivatives
(
const scalar time,
const scalarField& c,
const label li,
scalarField& dcdt
) const
{
const bool reduced = mechRed_->active();
const scalar T = c[this->nSpecie_];
const scalar p = c[this->nSpecie_ + 1];
if (reduced)
{
// When using DAC, the ODE solver submit a reduced set of species the
// complete set is used and only the species in the simplified mechanism
// are updated
this->c_ = completeC_;
// Update the concentration of the species in the simplified mechanism
// the other species remain the same and are used only for third-body
// efficiencies
for (label i=0; ic_[simplifiedToCompleteIndex_[i]] = max(c[i], 0);
}
}
else
{
for (label i=0; inSpecie(); i++)
{
this->c_[i] = max(c[i], 0);
}
}
omega(p, T, this->c_, li, dcdt);
// Constant pressure
// dT/dt = ...
scalar rho = 0;
for (label i=0; ic_.size(); i++)
{
const scalar W = this->specieThermo_[i].W();
rho += W*this->c_[i];
}
scalar cp = 0;
for (label i=0; ic_.size(); i++)
{
// cp function returns [J/kmol/K]
cp += this->c_[i]*this->specieThermo_[i].cp(p, T);
}
cp /= rho;
// When mechanism reduction is active
// dT is computed on the reduced set since dcdt is null
// for species not involved in the simplified mechanism
scalar dT = 0;
for (label i=0; inSpecie_; i++)
{
label si;
if (reduced)
{
si = simplifiedToCompleteIndex_[i];
}
else
{
si = i;
}
// ha function returns [J/kmol]
const scalar hi = this->specieThermo_[si].ha(p, T);
dT += hi*dcdt[i];
}
dT /= rho*cp;
dcdt[this->nSpecie_] = -dT;
// dp/dt = ...
dcdt[this->nSpecie_ + 1] = 0;
}
template
void Foam::TDACChemistryModel::jacobian
(
const scalar t,
const scalarField& c,
const label li,
scalarField& dcdt,
scalarSquareMatrix& J
) const
{
const bool reduced = mechRed_->active();
// If the mechanism reduction is active, the computed Jacobian
// is compact (size of the reduced set of species)
// but according to the information of the complete set
// (i.e. for the third-body efficiencies)
const scalar T = c[this->nSpecie_];
const scalar p = c[this->nSpecie_ + 1];
if (reduced)
{
this->c_ = completeC_;
for (label i=0; ic_[simplifiedToCompleteIndex_[i]] = max(c[i], 0);
}
}
else
{
forAll(this->c_, i)
{
this->c_[i] = max(c[i], 0);
}
}
J = Zero;
dcdt = Zero;
scalarField hi(this->c_.size());
scalarField cpi(this->c_.size());
forAll(hi, i)
{
hi[i] = this->specieThermo_[i].ha(p, T);
cpi[i] = this->specieThermo_[i].cp(p, T);
}
scalar omegaI = 0;
forAll(this->reactions_, ri)
{
if (!reactionsDisabled_[ri])
{
const Reaction& R = this->reactions_[ri];
scalar kfwd, kbwd;
R.dwdc
(
p,
T,
this->c_,
li,
J,
dcdt,
omegaI,
kfwd,
kbwd,
reduced,
completeToSimplifiedIndex_
);
R.dwdT
(
p,
T,
this->c_,
li,
omegaI,
kfwd,
kbwd,
J,
reduced,
completeToSimplifiedIndex_,
this->nSpecie_
);
}
}
// The species derivatives of the temperature term are partially computed
// while computing dwdc, they are completed hereunder:
scalar cpMean = 0;
scalar dcpdTMean = 0;
forAll(this->c_, i)
{
cpMean += this->c_[i]*cpi[i]; // J/(m^3 K)
// Already multiplied by rho
dcpdTMean += this->c_[i]*this->specieThermo_[i].dcpdT(p, T);
}
scalar dTdt = 0;
forAll(hi, i)
{
if (reduced)
{
const label si = completeToSimplifiedIndex_[i];
if (si != -1)
{
dTdt += hi[i]*dcdt[si]; // J/(m^3 s)
}
}
else
{
dTdt += hi[i]*dcdt[i]; // J/(m^3 s)
}
}
dTdt /= -cpMean; // K/s
dcdt[this->nSpecie_] = dTdt;
for (label i = 0; i < this->nSpecie_; i++)
{
J(this->nSpecie_, i) = 0;
for (label j = 0; j < this->nSpecie_; j++)
{
const label sj = reduced ? simplifiedToCompleteIndex_[j] : j;
J(this->nSpecie_, i) += hi[sj]*J(j, i);
}
const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
J(this->nSpecie_, i) += cpi[si]*dTdt; // J/(mol s)
J(this->nSpecie_, i) /= -cpMean; // K/s / (mol/m^3)
}
// ddT of dTdt
J(this->nSpecie_, this->nSpecie_) = 0;
for (label i = 0; i < this->nSpecie_; i++)
{
const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
J(this->nSpecie_, this->nSpecie_) +=
cpi[si]*dcdt[i]
+ hi[si]*J(i, this->nSpecie_);
}
J(this->nSpecie_, this->nSpecie_) += dTdt*dcpdTMean;
J(this->nSpecie_, this->nSpecie_) /= -cpMean;
J(this->nSpecie_, this->nSpecie_) += dTdt/T;
}
template
template
Foam::scalar Foam::TDACChemistryModel::solve
(
const DeltaTType& deltaT
)
{
// Increment counter of time-step
timeSteps_++;
const bool reduced = mechRed_->active();
label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
const basicSpecieMixture& composition = this->thermo().composition();
// CPU time analysis
const clockTime clockTime_= clockTime();
clockTime_.timeIncrement();
scalar reduceMechCpuTime_ = 0;
scalar addNewLeafCpuTime_ = 0;
scalar growCpuTime_ = 0;
scalar solveChemistryCpuTime_ = 0;
scalar searchISATCpuTime_ = 0;
this->resetTabulationResults();
// Average number of active species
scalar nActiveSpecies = 0;
scalar nAvg = 0;
BasicChemistryModel::correct();
scalar deltaTMin = great;
if (!this->chemistry_)
{
return deltaTMin;
}
const volScalarField rho
(
IOobject
(
"rho",
this->time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->thermo().rho()
);
const scalarField& T = this->thermo().T();
const scalarField& p = this->thermo().p();
scalarField c(this->nSpecie_);
scalarField c0(this->nSpecie_);
// Composition vector (Yi, T, p)
scalarField phiq(this->nEqns() + nAdditionalEqn);
scalarField Rphiq(this->nEqns() + nAdditionalEqn);
forAll(rho, celli)
{
const scalar rhoi = rho[celli];
scalar pi = p[celli];
scalar Ti = T[celli];
for (label i=0; inSpecie_; i++)
{
c[i] = rhoi*this->Y_[i][celli]/this->specieThermo_[i].W();
c0[i] = c[i];
phiq[i] = this->Y()[i][celli];
}
phiq[this->nSpecie()]=Ti;
phiq[this->nSpecie() + 1]=pi;
if (tabulation_->variableTimeStep())
{
phiq[this->nSpecie() + 2] = deltaT[celli];
}
// Initialise time progress
scalar timeLeft = deltaT[celli];
// Not sure if this is necessary
Rphiq = Zero;
clockTime_.timeIncrement();
// When tabulation is active (short-circuit evaluation for retrieve)
// It first tries to retrieve the solution of the system with the
// information stored through the tabulation method
if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
{
// Retrieved solution stored in Rphiq
for (label i=0; inSpecie(); i++)
{
c[i] = rhoi*Rphiq[i]/this->specieThermo_[i].W();
}
searchISATCpuTime_ += clockTime_.timeIncrement();
}
// This position is reached when tabulation is not used OR
// if the solution is not retrieved.
// In the latter case, it adds the information to the tabulation
// (it will either expand the current data or add a new stored point).
else
{
searchISATCpuTime_ += clockTime_.timeIncrement();
if (reduced)
{
// Reduce mechanism change the number of species (only active)
mechRed_->reduceMechanism(pi, Ti, c, celli);
nActiveSpecies += mechRed_->NsSimp();
nAvg++;
reduceMechCpuTime_ += timeIncr;
}
// Calculate the chemical source terms
while (timeLeft > small)
{
scalar dt = timeLeft;
if (reduced)
{
// completeC_ used in the overridden ODE methods
// to update only the active species
completeC_ = c;
// Solve the reduced set of ODE
this->solve
(
pi,
Ti,
simplifiedC_,
celli,
dt,
this->deltaTChem_[celli]
);
for (label i=0; isolve(pi, Ti, c, celli, dt, this->deltaTChem_[celli]);
}
timeLeft -= dt;
}
{
solveChemistryCpuTime_ += timeIncr;
}
// If tabulation is used, we add the information computed here to
// the stored points (either expand or add)
if (tabulation_->active())
{
forAll(c, i)
{
Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W();
}
if (tabulation_->variableTimeStep())
{
Rphiq[Rphiq.size()-3] = Ti;
Rphiq[Rphiq.size()-2] = pi;
Rphiq[Rphiq.size()-1] = deltaT[celli];
}
else
{
Rphiq[Rphiq.size()-2] = Ti;
Rphiq[Rphiq.size()-1] = pi;
}
label growOrAdd =
tabulation_->add(phiq, Rphiq, celli, rhoi, deltaT[celli]);
if (growOrAdd)
{
this->setTabulationResultsAdd(celli);
addNewLeafCpuTime_ += clockTime_.timeIncrement();
}
else
{
this->setTabulationResultsGrow(celli);
growCpuTime_ += clockTime_.timeIncrement();
}
}
// When operations are done and if mechanism reduction is active,
// the number of species (which also affects nEqns) is set back
// to the total number of species (stored in the mechRed object)
if (reduced)
{
this->nSpecie_ = mechRed_->nSpecie();
}
deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
this->deltaTChem_[celli] =
min(this->deltaTChem_[celli], this->deltaTChemMax_);
}
// Set the RR vector (used in the solver)
for (label i=0; inSpecie_; i++)
{
this->RR_[i][celli] =
(c[i] - c0[i])*this->specieThermo_[i].W()/deltaT[celli];
}
}
if (mechRed_->log() || tabulation_->log())
{
cpuSolveFile_()
<< this->time().timeOutputValue()
<< " " << solveChemistryCpuTime_ << endl;
}
if (mechRed_->log())
{
cpuReduceFile_()
<< this->time().timeOutputValue()
<< " " << reduceMechCpuTime_ << endl;
}
if (tabulation_->active())
{
// Every time-step, look if the tabulation should be updated
tabulation_->update();
// Write the performance of the tabulation
tabulation_->writePerformance();
if (tabulation_->log())
{
cpuRetrieveFile_()
<< this->time().timeOutputValue()
<< " " << searchISATCpuTime_ << endl;
cpuGrowFile_()
<< this->time().timeOutputValue()
<< " " << growCpuTime_ << endl;
cpuAddFile_()
<< this->time().timeOutputValue()
<< " " << addNewLeafCpuTime_ << endl;
}
}
if (reduced && nAvg && mechRed_->log())
{
// Write average number of species
nActiveSpeciesFile_()
<< this->time().timeOutputValue()
<< " " << nActiveSpecies/nAvg << endl;
}
if (Pstream::parRun())
{
List active(composition.active());
Pstream::listCombineGather(active, orEqOp());
Pstream::listCombineScatter(active);
forAll(active, i)
{
if (active[i])
{
composition.setActive(i);
}
}
}
return deltaTMin;
}
template
Foam::scalar Foam::TDACChemistryModel::solve
(
const scalar deltaT
)
{
// Don't allow the time-step to change more than a factor of 2
return min
(
this->solve>(UniformField(deltaT)),
2*deltaT
);
}
template
Foam::scalar Foam::TDACChemistryModel::solve
(
const scalarField& deltaT
)
{
return this->solve(deltaT);
}
template
void Foam::TDACChemistryModel::
setTabulationResultsAdd
(
const label celli
)
{
tabulationResults_[celli] = 0;
}
template
void Foam::TDACChemistryModel::
setTabulationResultsGrow(const label celli)
{
tabulationResults_[celli] = 1;
}
template
void Foam::TDACChemistryModel::
setTabulationResultsRetrieve
(
const label celli
)
{
tabulationResults_[celli] = 2;
}
// ************************************************************************* //