Commit 7e22440d authored by Henry Weller's avatar Henry Weller
Browse files

TDACChemistryModel: Added support for variable time-step and LTS in ISAT

New reactingFoam tutorial counterFlowFlame2DLTS_GRI_TDAC demonstrates this new
functionality.

Additionally the ISAT table growth algorithm has been further optimized
providing an overall speedup of between 15% and 38% for the tests run so far.

Updates to TDAC and ISAT provided by Francesco Contino.

Implementation updated and integrated into OpenFOAM-dev by
Henry G. Weller, CFD Direct Ltd with the help of Francesco Contino.

Original code providing all algorithms for chemistry reduction and
tabulation contributed by Francesco Contino, Tommaso Lucchini, Gianluca
D’Errico, Hervé Jeanmart, Nicolas Bourgeois and Stéphane Backaert.
parent 126125c1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -37,12 +37,26 @@ Foam::TDACChemistryModel<CompType, ThermoType>::TDACChemistryModel
)
:
chemistryModel<CompType, ThermoType>(mesh, phaseName),
timeSteps_(0),
NsDAC_(this->nSpecie_),
completeC_(this->nSpecie_,0.0),
completeC_(this->nSpecie_, 0),
reactionsDisabled_(this->reactions_.size(), false),
completeToSimplifiedIndex_(this->nSpecie_,-1),
specieComp_(this->nSpecie_),
completeToSimplifiedIndex_(this->nSpecie_, -1),
simplifiedToCompleteIndex_(this->nSpecie_),
specieComp_(this->nSpecie_)
tabulationResults_
(
IOobject
(
"TabulationResults",
this->time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
0
)
{
basicMultiComponentMixture& composition = this->thermo().composition();
......@@ -53,7 +67,7 @@ Foam::TDACChemistryModel<CompType, ThermoType>::TDACChemistryModel
dynamicCast<const reactingMixture<ThermoType>&>(this->thermo())
.specieComposition();
forAll(specieComp_,i)
forAll(specieComp_, i)
{
specieComp_[i] = specComp[this->Y()[i].name()];
}
......@@ -579,8 +593,13 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
const DeltaTType& deltaT
)
{
// Increment counter of time-step
timeSteps_++;
const bool reduced = mechRed_->active();
label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1:0);
basicMultiComponentMixture& composition = this->thermo().composition();
// CPU time analysis
......@@ -625,9 +644,9 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
scalarField c0(this->nSpecie_);
// Composition vector (Yi, T, p)
scalarField phiq(this->nEqns());
scalarField phiq(this->nEqns() + nAdditionalEqn);
scalarField Rphiq(this->nEqns());
scalarField Rphiq(this->nEqns() + nAdditionalEqn);
forAll(rho, celli)
{
......@@ -643,6 +662,11 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
}
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];
......@@ -668,7 +692,7 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
// 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 poin).
// (it will either expand the current data or add a new stored point).
else
{
clockTime_.timeIncrement();
......@@ -720,12 +744,28 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
{
Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W();
}
Rphiq[Rphiq.size()-2] = Ti;
Rphiq[Rphiq.size()-1] = pi;
tabulation_->add(phiq, Rphiq, rhoi);
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, rhoi, deltaT[celli]);
if (growOrAdd)
{
this->setTabulationResultsAdd(celli);
}
else
{
this->setTabulationResultsGrow(celli);
}
}
addNewLeafCpuTime_ += clockTime_.timeIncrement();
// When operations are done and if mechanism reduction is active,
......@@ -840,4 +880,35 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
}
template<class CompType, class ThermoType>
void Foam::TDACChemistryModel<CompType, ThermoType>::setTabulationResultsAdd
(
const label celli
)
{
tabulationResults_[celli] = 0.0;
}
template<class CompType, class ThermoType>
void Foam::TDACChemistryModel<CompType, ThermoType>::setTabulationResultsGrow
(
const label celli
)
{
tabulationResults_[celli] = 1.0;
}
template<class CompType, class ThermoType>
void Foam::TDACChemistryModel<CompType, ThermoType>::
setTabulationResultsRetrieve
(
const label celli
)
{
tabulationResults_[celli] = 2.0;
}
// ************************************************************************* //
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -85,14 +85,16 @@ class TDACChemistryModel
{
// Private member data
label timeSteps_;
// Mechanism reduction
label NsDAC_;
scalarField completeC_;
scalarField simplifiedC_;
Field<bool> reactionsDisabled_;
List<List<specieElement>> specieComp_;
Field<label> completeToSimplifiedIndex_;
DynamicList<label> simplifiedToCompleteIndex_;
List<List<specieElement>> specieComp_;
autoPtr<chemistryReductionMethod<CompType, ThermoType>> mechRed_;
// Tabulation
......@@ -113,6 +115,12 @@ class TDACChemistryModel
//- Log file for average time spent solving the chemistry
autoPtr<OFstream> cpuSolveFile_;
// Field containing information about tabulation:
// 0 -> add (direct integration)
// 1 -> grow
// 2 -> retrieve
volScalarField tabulationResults_;
// Private Member Functions
......@@ -151,6 +159,12 @@ public:
// Member Functions
inline label timeSteps() const
{
return timeSteps_;
}
//- Create and return a TDAC log file of the given name
inline autoPtr<OFstream> logFile(const word& name) const;
......@@ -256,6 +270,17 @@ public:
inline autoPtr<chemistryReductionMethod<CompType, ThermoType>>&
mechRed();
tmp<volScalarField> tabulationResults() const
{
return tabulationResults_;
}
void setTabulationResultsAdd(const label celli);
void setTabulationResultsGrow(const label celli);
void setTabulationResultsRetrieve(const label celli);
};
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -41,16 +41,11 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
chemistry
),
chemisTree_(chemistry, this->coeffsDict_),
scaleFactor_(chemistry.nEqns(), 1.0),
scaleFactor_(chemistry.nEqns() + ((this->variableTimeStep()) ? 1 : 0), 1),
runTime_(chemistry.time()),
chPMaxLifeTime_
(
this->coeffsDict_.lookupOrDefault
(
"chPMaxLifeTime",
(runTime_.endTime().value()-runTime_.startTime().value())
/runTime_.deltaT().value()
)
this->coeffsDict_.lookupOrDefault("chPMaxLifeTime", INT_MAX)
),
maxGrowth_(this->coeffsDict_.lookupOrDefault("maxGrowth", INT_MAX)),
checkEntireTreeInterval_
......@@ -104,6 +99,36 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
}
scaleFactor_[Ysize] = readScalar(scaleDict.lookup("Temperature"));
scaleFactor_[Ysize+1] = readScalar(scaleDict.lookup("Pressure"));
if (this->variableTimeStep())
{
scaleFactor_[Ysize + 2] = readScalar(scaleDict.lookup("deltaT"));
}
else
{
// When the variableTimeStep option is false, if the application
// has variable time step activated, the maximum lifetime of a
// chemPoint should be 1 time step.
bool adjustTimeStep =
runTime_.controlDict().lookupOrDefault("adjustTimeStep", false);
if (chPMaxLifeTime_ > 1 && adjustTimeStep)
{
WarningInFunction
<< " variableTimeStep is not activate for ISAT while"
<< " the time step might be adjusted by the application."
<< nl
<< " This might lead to errors in the chemistry." << nl
<< " To avoid this warning either set chPMaxLifeTime to 1"
<< " or activate variableTimeStep." << endl;
}
}
}
if (this->variableTimeStep())
{
nAdditionalEqns_ = 3;
}
else
{
nAdditionalEqns_ = 2;
}
if (this->log())
......@@ -190,16 +215,16 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
scalarField& Rphiq
)
{
label nEqns = this->chemistry_.nEqns(); // Full set of species
label nEqns = this->chemistry_.nEqns(); // Species, T, p
bool mechRedActive = this->chemistry_.mechRed()->active();
Rphiq = phi0->Rphi();
scalarField dphi(phiq-phi0->phi());
const scalarSquareMatrix& gradientsMatrix = phi0->A();
List<label>& completeToSimplified(phi0->completeToSimplifiedIndex());
// Rphiq[i] = Rphi0[i]+A(i, j)dphi[j]
// Rphiq[i]=Rphi0[i]+A(i, j)dphi[j]
// where Aij is dRi/dphi_j
for (label i=0; i<nEqns-2; i++)
for (label i=0; i<nEqns-nAdditionalEqns_; i++)
{
if (mechRedActive)
{
......@@ -216,9 +241,18 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
}
}
Rphiq[i] +=
gradientsMatrix(si, phi0->nActiveSpecies())*dphi[nEqns-2];
gradientsMatrix(si, phi0->nActiveSpecies())*dphi[nEqns - 2];
Rphiq[i] +=
gradientsMatrix(si, phi0->nActiveSpecies()+1)*dphi[nEqns-1];
gradientsMatrix(si, phi0->nActiveSpecies() + 1)
*dphi[nEqns - 1];
if (this->variableTimeStep())
{
Rphiq[i] +=
gradientsMatrix(si, phi0->nActiveSpecies() + 2)
*dphi[nEqns];
}
// As we use an approximation of A, Rphiq should be checked for
// negative values
Rphiq[i] = max(0.0,Rphiq[i]);
......@@ -260,10 +294,11 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::grow
// Raise a flag when the chemPoint used has been grown more than the
// allowed number of time
if (!phi0->toRemove() && phi0->nGrowth() > maxGrowth_)
if (phi0->nGrowth() > maxGrowth_)
{
cleaningRequired_ = true;
phi0->toRemove() = true;
return false;
}
// If the solution RphiQ is still within the tolerance we try to grow it
......@@ -294,14 +329,10 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
{
chemPointISAT<CompType, ThermoType>* xtmp =
chemisTree_.treeSuccessor(x);
// timeOutputValue returns timeToUserTime(value()), therefore, it should
// be compare with timeToUserTime(deltaT)
scalar elapsedTime = runTime_.timeOutputValue() - x->timeTag();
scalar maxElapsedTime =
chPMaxLifeTime_
* runTime_.timeToUserTime(runTime_.deltaTValue());
if ((elapsedTime > maxElapsedTime) || (x->nGrowth() > maxGrowth_))
scalar elapsedTimeSteps = this->chemistry_.timeSteps() - x->timeTag();
if ((elapsedTimeSteps > chPMaxLifeTime_) || (x->nGrowth() > maxGrowth_))
{
chemisTree_.deleteLeaf(x);
treeModified = true;
......@@ -334,13 +365,13 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
(
scalarSquareMatrix& A,
const scalarField& Rphiq,
const scalar rhoi
const scalar rhoi,
const scalar dt
)
{
scalar dt = runTime_.deltaTValue();
bool mechRedActive = this->chemistry_.mechRed()->active();
label speciesNumber = this->chemistry_.nSpecie();
scalarField Rcq(this->chemistry_.nEqns());
scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
for (label i=0; i<speciesNumber; i++)
{
label s2c = i;
......@@ -350,8 +381,12 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
}
Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
}
Rcq[speciesNumber] = Rphiq[Rphiq.size()-2];
Rcq[speciesNumber+1] = Rphiq[Rphiq.size()-1];
Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
if (this->variableTimeStep())
{
Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
}
// Aaa is computed implicitely,
// A is given by A = C(psi0, t0+dt), where C is obtained through solving
......@@ -399,6 +434,10 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
// For temperature and pressure, only unity on the diagonal
A(speciesNumber, speciesNumber) = 1;
A(speciesNumber+1, speciesNumber+1) = 1;
if (this->variableTimeStep())
{
A[speciesNumber + 2][speciesNumber + 2] = 1;
}
// Inverse of (I-dt*J(psi(t0+dt)))
LUscalarMatrix LUA(A);
......@@ -463,20 +502,18 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
if (retrieved)
{
scalar elapsedTime =
runTime_.timeOutputValue() - phi0->timeTag();
scalar maxElapsedTime =
chPMaxLifeTime_
* runTime_.timeToUserTime(runTime_.deltaTValue());
phi0->increaseNumRetrieve();
scalar elapsedTimeSteps =
this->chemistry_.timeSteps() - phi0->timeTag();
// Raise a flag when the chemPoint has been used more than the allowed
// number of time steps
if (elapsedTime > maxElapsedTime && !phi0->toRemove())
if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->toRemove())
{
cleaningRequired_ = true;
phi0->toRemove() = true;
}
lastSearch_->lastTimeUsed() = runTime_.timeOutputValue();
lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
addToMRU(phi0);
calcNewC(phi0,phiq, Rphiq);
nRetrieved_++;
......@@ -492,13 +529,15 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
template<class CompType, class ThermoType>
bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
(
const scalarField& phiq,
const scalarField& Rphiq,
const scalar rho
const scalar rho,
const scalar deltaT
)
{
label growthOrAddFlag = 1;
// If lastSearch_ holds a valid pointer to a chemPoint AND the growPoints_
// option is on, the code first tries to grow the point hold by lastSearch_
if (lastSearch_ && growPoints_)
......@@ -506,13 +545,12 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
if (grow(lastSearch_,phiq, Rphiq))
{
nGrowth_++;
// the structure of the tree is not modified, return false
return false;
growthOrAddFlag = 0;
//the structure of the tree is not modified, return false
return growthOrAddFlag;
}
}
bool treeCleanedOrCleared(false);
// If the code reach this point, it is either because lastSearch_ is not
// valid, OR because growPoints_ is not on, OR because the grow operation
// has failed. In the three cases, a new point is added to the tree.
......@@ -567,16 +605,12 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
// The structure has been changed, it will force the binary tree to
// perform a new search and find the most appropriate point still stored
lastSearch_ = nullptr;
// Either cleanAndBalance has changed the tree or it has been cleared
// in any case treeCleanedOrCleared should be set to true
treeCleanedOrCleared = true;
}
// Compute the A matrix needed to store the chemPoint.
label ASize = this->chemistry_.nEqns(); // Reduced when mechRed is active
label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
scalarSquareMatrix A(ASize, Zero);
computeA(A, Rphiq, rho);
computeA(A, Rphiq, rho, deltaT);
chemisTree().insertNewLeaf
(
......@@ -591,7 +625,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
nAdd_++;
return treeCleanedOrCleared;
return growthOrAddFlag;
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -112,6 +112,9 @@ class ISAT
bool cleaningRequired_;
//- Number of equations in addition to the species eqs.
label nAdditionalEqns_;
// Private Member Functions
......@@ -168,7 +171,8 @@ class ISAT
(
scalarSquareMatrix& A,
const scalarField& Rphiq,
const scalar rho
const scalar rho,
const scalar dt
);
......@@ -224,11 +228,12 @@ public:
// This function can grow an existing point or add a new leaf to the
// binary tree Input : phiq the new composition to store Rphiq the
// mapping of the new composition point
virtual bool add
virtual label add
(
const scalarField& phiq,
const scalarField& Rphiq,
const scalar rho
const scalar rho,
const scalar deltaT
);
virtual bool update()
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -36,8 +36,11 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
leafRight_(nullptr),
nodeLeft_(nullptr),
nodeRight_(nullptr),
parent_(nullptr)
{}
parent_(nullptr),
variableTimeStep_(false),
nAdditionalEqns_(0)
{
}
template<class CompType, class ThermoType>
......@@ -53,8 +56,18 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
nodeLeft_(nullptr),
nodeRight_(nullptr),
parent_(parent),
v_(elementLeft->completeSpaceSize(),0.0)