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.
This commit is contained in:
Henry Weller
2017-01-07 16:29:15 +00:00
parent 78a396430b
commit 923350fa6e
36 changed files with 8750 additions and 205 deletions

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -37,12 +37,26 @@ Foam::TDACChemistryModel<CompType, ThermoType>::TDACChemistryModel
) )
: :
chemistryModel<CompType, ThermoType>(mesh, phaseName), chemistryModel<CompType, ThermoType>(mesh, phaseName),
timeSteps_(0),
NsDAC_(this->nSpecie_), NsDAC_(this->nSpecie_),
completeC_(this->nSpecie_,0.0), completeC_(this->nSpecie_, 0),
reactionsDisabled_(this->reactions_.size(), false), reactionsDisabled_(this->reactions_.size(), false),
completeToSimplifiedIndex_(this->nSpecie_,-1), specieComp_(this->nSpecie_),
completeToSimplifiedIndex_(this->nSpecie_, -1),
simplifiedToCompleteIndex_(this->nSpecie_), 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(); basicMultiComponentMixture& composition = this->thermo().composition();
@ -53,7 +67,7 @@ Foam::TDACChemistryModel<CompType, ThermoType>::TDACChemistryModel
dynamicCast<const reactingMixture<ThermoType>&>(this->thermo()) dynamicCast<const reactingMixture<ThermoType>&>(this->thermo())
.specieComposition(); .specieComposition();
forAll(specieComp_,i) forAll(specieComp_, i)
{ {
specieComp_[i] = specComp[this->Y()[i].name()]; specieComp_[i] = specComp[this->Y()[i].name()];
} }
@ -579,8 +593,13 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
const DeltaTType& deltaT const DeltaTType& deltaT
) )
{ {
// Increment counter of time-step
timeSteps_++;
const bool reduced = mechRed_->active(); const bool reduced = mechRed_->active();
label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1:0);
basicMultiComponentMixture& composition = this->thermo().composition(); basicMultiComponentMixture& composition = this->thermo().composition();
// CPU time analysis // CPU time analysis
@ -625,9 +644,9 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
scalarField c0(this->nSpecie_); scalarField c0(this->nSpecie_);
// Composition vector (Yi, T, p) // 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) forAll(rho, celli)
{ {
@ -643,6 +662,11 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
} }
phiq[this->nSpecie()]=Ti; phiq[this->nSpecie()]=Ti;
phiq[this->nSpecie()+1]=pi; phiq[this->nSpecie()+1]=pi;
if (tabulation_->variableTimeStep())
{
phiq[this->nSpecie()+2] = deltaT[celli];
}
// Initialise time progress // Initialise time progress
scalar timeLeft = deltaT[celli]; 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 // This position is reached when tabulation is not used OR
// if the solution is not retrieved. // if the solution is not retrieved.
// In the latter case, it adds the information to the tabulation // 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 else
{ {
clockTime_.timeIncrement(); clockTime_.timeIncrement();
@ -720,12 +744,28 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
{ {
Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W(); Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W();
} }
if (tabulation_->variableTimeStep())
Rphiq[Rphiq.size()-2] = Ti; {
Rphiq[Rphiq.size()-1] = pi; Rphiq[Rphiq.size()-3] = Ti;
tabulation_->add(phiq, Rphiq, rhoi); 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(); addNewLeafCpuTime_ += clockTime_.timeIncrement();
// When operations are done and if mechanism reduction is active, // 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;
}
// ************************************************************************* // // ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -85,14 +85,16 @@ class TDACChemistryModel
{ {
// Private member data // Private member data
label timeSteps_;
// Mechanism reduction // Mechanism reduction
label NsDAC_; label NsDAC_;
scalarField completeC_; scalarField completeC_;
scalarField simplifiedC_; scalarField simplifiedC_;
Field<bool> reactionsDisabled_; Field<bool> reactionsDisabled_;
List<List<specieElement>> specieComp_;
Field<label> completeToSimplifiedIndex_; Field<label> completeToSimplifiedIndex_;
DynamicList<label> simplifiedToCompleteIndex_; DynamicList<label> simplifiedToCompleteIndex_;
List<List<specieElement>> specieComp_;
autoPtr<chemistryReductionMethod<CompType, ThermoType>> mechRed_; autoPtr<chemistryReductionMethod<CompType, ThermoType>> mechRed_;
// Tabulation // Tabulation
@ -113,6 +115,12 @@ class TDACChemistryModel
//- Log file for average time spent solving the chemistry //- Log file for average time spent solving the chemistry
autoPtr<OFstream> cpuSolveFile_; autoPtr<OFstream> cpuSolveFile_;
// Field containing information about tabulation:
// 0 -> add (direct integration)
// 1 -> grow
// 2 -> retrieve
volScalarField tabulationResults_;
// Private Member Functions // Private Member Functions
@ -151,6 +159,12 @@ public:
// Member Functions // Member Functions
inline label timeSteps() const
{
return timeSteps_;
}
//- Create and return a TDAC log file of the given name //- Create and return a TDAC log file of the given name
inline autoPtr<OFstream> logFile(const word& name) const; inline autoPtr<OFstream> logFile(const word& name) const;
@ -256,6 +270,17 @@ public:
inline autoPtr<chemistryReductionMethod<CompType, ThermoType>>& inline autoPtr<chemistryReductionMethod<CompType, ThermoType>>&
mechRed(); mechRed();
tmp<volScalarField> tabulationResults() const
{
return tabulationResults_;
}
void setTabulationResultsAdd(const label celli);
void setTabulationResultsGrow(const label celli);
void setTabulationResultsRetrieve(const label celli);
}; };

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -41,16 +41,11 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
chemistry chemistry
), ),
chemisTree_(chemistry, this->coeffsDict_), chemisTree_(chemistry, this->coeffsDict_),
scaleFactor_(chemistry.nEqns(), 1.0), scaleFactor_(chemistry.nEqns() + ((this->variableTimeStep()) ? 1 : 0), 1),
runTime_(chemistry.time()), runTime_(chemistry.time()),
chPMaxLifeTime_ chPMaxLifeTime_
( (
this->coeffsDict_.lookupOrDefault this->coeffsDict_.lookupOrDefault("chPMaxLifeTime", INT_MAX)
(
"chPMaxLifeTime",
(runTime_.endTime().value()-runTime_.startTime().value())
/runTime_.deltaT().value()
)
), ),
maxGrowth_(this->coeffsDict_.lookupOrDefault("maxGrowth", INT_MAX)), maxGrowth_(this->coeffsDict_.lookupOrDefault("maxGrowth", INT_MAX)),
checkEntireTreeInterval_ checkEntireTreeInterval_
@ -104,6 +99,36 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
} }
scaleFactor_[Ysize] = readScalar(scaleDict.lookup("Temperature")); scaleFactor_[Ysize] = readScalar(scaleDict.lookup("Temperature"));
scaleFactor_[Ysize+1] = readScalar(scaleDict.lookup("Pressure")); 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()) if (this->log())
@ -190,16 +215,16 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
scalarField& Rphiq 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(); bool mechRedActive = this->chemistry_.mechRed()->active();
Rphiq = phi0->Rphi(); Rphiq = phi0->Rphi();
scalarField dphi(phiq-phi0->phi()); scalarField dphi(phiq-phi0->phi());
const scalarSquareMatrix& gradientsMatrix = phi0->A(); const scalarSquareMatrix& gradientsMatrix = phi0->A();
List<label>& completeToSimplified(phi0->completeToSimplifiedIndex()); 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 // where Aij is dRi/dphi_j
for (label i=0; i<nEqns-2; i++) for (label i=0; i<nEqns-nAdditionalEqns_; i++)
{ {
if (mechRedActive) if (mechRedActive)
{ {
@ -216,9 +241,18 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
} }
} }
Rphiq[i] += Rphiq[i] +=
gradientsMatrix(si, phi0->nActiveSpecies())*dphi[nEqns-2]; gradientsMatrix(si, phi0->nActiveSpecies())*dphi[nEqns - 2];
Rphiq[i] += 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 // As we use an approximation of A, Rphiq should be checked for
// negative values // negative values
Rphiq[i] = max(0.0,Rphiq[i]); 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 // Raise a flag when the chemPoint used has been grown more than the
// allowed number of time // allowed number of time
if (!phi0->toRemove() && phi0->nGrowth() > maxGrowth_) if (phi0->nGrowth() > maxGrowth_)
{ {
cleaningRequired_ = true; cleaningRequired_ = true;
phi0->toRemove() = true; phi0->toRemove() = true;
return false;
} }
// If the solution RphiQ is still within the tolerance we try to grow it // 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 = chemPointISAT<CompType, ThermoType>* xtmp =
chemisTree_.treeSuccessor(x); 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); chemisTree_.deleteLeaf(x);
treeModified = true; treeModified = true;
@ -334,13 +365,13 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
( (
scalarSquareMatrix& A, scalarSquareMatrix& A,
const scalarField& Rphiq, const scalarField& Rphiq,
const scalar rhoi const scalar rhoi,
const scalar dt
) )
{ {
scalar dt = runTime_.deltaTValue();
bool mechRedActive = this->chemistry_.mechRed()->active(); bool mechRedActive = this->chemistry_.mechRed()->active();
label speciesNumber = this->chemistry_.nSpecie(); label speciesNumber = this->chemistry_.nSpecie();
scalarField Rcq(this->chemistry_.nEqns()); scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
for (label i=0; i<speciesNumber; i++) for (label i=0; i<speciesNumber; i++)
{ {
label s2c = 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[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
} }
Rcq[speciesNumber] = Rphiq[Rphiq.size()-2]; Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
Rcq[speciesNumber+1] = Rphiq[Rphiq.size()-1]; Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
if (this->variableTimeStep())
{
Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
}
// Aaa is computed implicitely, // Aaa is computed implicitely,
// A is given by A = C(psi0, t0+dt), where C is obtained through solving // 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 // For temperature and pressure, only unity on the diagonal
A(speciesNumber, speciesNumber) = 1; A(speciesNumber, speciesNumber) = 1;
A(speciesNumber+1, speciesNumber+1) = 1; A(speciesNumber+1, speciesNumber+1) = 1;
if (this->variableTimeStep())
{
A[speciesNumber + 2][speciesNumber + 2] = 1;
}
// Inverse of (I-dt*J(psi(t0+dt))) // Inverse of (I-dt*J(psi(t0+dt)))
LUscalarMatrix LUA(A); LUscalarMatrix LUA(A);
@ -463,20 +502,18 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
if (retrieved) if (retrieved)
{ {
scalar elapsedTime = phi0->increaseNumRetrieve();
runTime_.timeOutputValue() - phi0->timeTag(); scalar elapsedTimeSteps =
scalar maxElapsedTime = this->chemistry_.timeSteps() - phi0->timeTag();
chPMaxLifeTime_
* runTime_.timeToUserTime(runTime_.deltaTValue());
// Raise a flag when the chemPoint has been used more than the allowed // Raise a flag when the chemPoint has been used more than the allowed
// number of time steps // number of time steps
if (elapsedTime > maxElapsedTime && !phi0->toRemove()) if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->toRemove())
{ {
cleaningRequired_ = true; cleaningRequired_ = true;
phi0->toRemove() = true; phi0->toRemove() = true;
} }
lastSearch_->lastTimeUsed() = runTime_.timeOutputValue(); lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
addToMRU(phi0); addToMRU(phi0);
calcNewC(phi0,phiq, Rphiq); calcNewC(phi0,phiq, Rphiq);
nRetrieved_++; nRetrieved_++;
@ -492,13 +529,15 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
template<class CompType, class ThermoType> 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& phiq,
const scalarField& Rphiq, 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_ // 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_ // option is on, the code first tries to grow the point hold by lastSearch_
if (lastSearch_ && growPoints_) if (lastSearch_ && growPoints_)
@ -506,13 +545,12 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
if (grow(lastSearch_,phiq, Rphiq)) if (grow(lastSearch_,phiq, Rphiq))
{ {
nGrowth_++; nGrowth_++;
// the structure of the tree is not modified, return false growthOrAddFlag = 0;
return false; //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 // 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 // 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. // 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 // The structure has been changed, it will force the binary tree to
// perform a new search and find the most appropriate point still stored // perform a new search and find the most appropriate point still stored
lastSearch_ = nullptr; 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. // 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); scalarSquareMatrix A(ASize, Zero);
computeA(A, Rphiq, rho); computeA(A, Rphiq, rho, deltaT);
chemisTree().insertNewLeaf chemisTree().insertNewLeaf
( (
@ -591,7 +625,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
nAdd_++; nAdd_++;
return treeCleanedOrCleared; return growthOrAddFlag;
} }

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -112,6 +112,9 @@ class ISAT
bool cleaningRequired_; bool cleaningRequired_;
//- Number of equations in addition to the species eqs.
label nAdditionalEqns_;
// Private Member Functions // Private Member Functions
@ -168,7 +171,8 @@ class ISAT
( (
scalarSquareMatrix& A, scalarSquareMatrix& A,
const scalarField& Rphiq, 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 // 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 // binary tree Input : phiq the new composition to store Rphiq the
// mapping of the new composition point // mapping of the new composition point
virtual bool add virtual label add
( (
const scalarField& phiq, const scalarField& phiq,
const scalarField& Rphiq, const scalarField& Rphiq,
const scalar rho const scalar rho,
const scalar deltaT
); );
virtual bool update() virtual bool update()

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -36,8 +36,11 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
leafRight_(nullptr), leafRight_(nullptr),
nodeLeft_(nullptr), nodeLeft_(nullptr),
nodeRight_(nullptr), nodeRight_(nullptr),
parent_(nullptr) parent_(nullptr),
{} variableTimeStep_(false),
nAdditionalEqns_(0)
{
}
template<class CompType, class ThermoType> template<class CompType, class ThermoType>
@ -53,8 +56,18 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
nodeLeft_(nullptr), nodeLeft_(nullptr),
nodeRight_(nullptr), nodeRight_(nullptr),
parent_(parent), parent_(parent),
v_(elementLeft->completeSpaceSize(),0.0) variableTimeStep_(elementLeft->variableTimeStep()),
v_(elementLeft->completeSpaceSize(), 0)
{ {
if (this->variableTimeStep_)
{
nAdditionalEqns_ = 3;
}
else
{
nAdditionalEqns_ = 2;
}
calcV(elementLeft, elementRight, v_); calcV(elementLeft, elementRight, v_);
a_ = calcA(elementLeft, elementRight); a_ = calcA(elementLeft, elementRight);
} }
@ -70,9 +83,22 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
nodeLeft_(bn->nodeLeft()), nodeLeft_(bn->nodeLeft()),
nodeRight_(bn->nodeRight()), nodeRight_(bn->nodeRight()),
parent_(bn->parent()), parent_(bn->parent()),
variableTimeStep_
(
this->coeffsDict_.lookupOrDefault("variableTimeStep", false)
),
v_(bn->v()), v_(bn->v()),
a_(bn->a()) a_(bn->a())
{} {
if (this->variableTimeStep_)
{
nAdditionalEqns_ = 3;
}
else
{
nAdditionalEqns_ = 2;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -102,7 +128,7 @@ Foam::binaryNode<CompType, ThermoType>::calcV
bool outOfIndexI = true; bool outOfIndexI = true;
if (mechReductionActive) if (mechReductionActive)
{ {
if (i<elementLeft->completeSpaceSize()-2) if (i<elementLeft->completeSpaceSize() - nAdditionalEqns_)
{ {
si = elementLeft->completeToSimplifiedIndex()[i]; si = elementLeft->completeToSimplifiedIndex()[i];
outOfIndexI = (si==-1); outOfIndexI = (si==-1);
@ -110,8 +136,9 @@ Foam::binaryNode<CompType, ThermoType>::calcV
else// temperature and pressure else// temperature and pressure
{ {
outOfIndexI = false; outOfIndexI = false;
label dif = i-(elementLeft->completeSpaceSize()-2); const label dif =
si = elementLeft->nActiveSpecies()+dif; i - (elementLeft->completeSpaceSize() - nAdditionalEqns_);
si = elementLeft->nActiveSpecies() + dif;
} }
} }
if (!mechReductionActive || (mechReductionActive && !(outOfIndexI))) if (!mechReductionActive || (mechReductionActive && !(outOfIndexI)))
@ -123,7 +150,7 @@ Foam::binaryNode<CompType, ThermoType>::calcV
bool outOfIndexJ = true; bool outOfIndexJ = true;
if (mechReductionActive) if (mechReductionActive)
{ {
if (j<elementLeft->completeSpaceSize()-2) if (j < elementLeft->completeSpaceSize() - nAdditionalEqns_)
{ {
sj = elementLeft->completeToSimplifiedIndex()[j]; sj = elementLeft->completeToSimplifiedIndex()[j];
outOfIndexJ = (sj==-1); outOfIndexJ = (sj==-1);
@ -131,8 +158,13 @@ Foam::binaryNode<CompType, ThermoType>::calcV
else else
{ {
outOfIndexJ = false; outOfIndexJ = false;
label dif = j-(elementLeft->completeSpaceSize()-2); const label dif =
sj = elementLeft->nActiveSpecies()+dif; j
- (
elementLeft->completeSpaceSize()
- nAdditionalEqns_
);
sj = elementLeft->nActiveSpecies() + dif;
} }
} }
if if

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -67,7 +67,14 @@ public:
//- Parent node //- Parent node
binaryNode<CompType, ThermoType>* parent_; binaryNode<CompType, ThermoType>* parent_;
//- Switch to allow variable time step (off by default)
bool variableTimeStep_;
//- Number of equations in addition to the species eqs.
label nAdditionalEqns_;
scalarField v_; scalarField v_;
scalar a_; scalar a_;
//- Compute vector v: //- Compute vector v:

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -66,7 +66,7 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
{ {
if ((n2ndSearch_ < max2ndSearch_) && (y!=nullptr)) if ((n2ndSearch_ < max2ndSearch_) && (y!=nullptr))
{ {
scalar vPhi=0.0; scalar vPhi = 0;
const scalarField& v = y->v(); const scalarField& v = y->v();
const scalar a = y->a(); const scalar a = y->a();
// compute v*phi // compute v*phi
@ -74,18 +74,18 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
{ {
vPhi += phiq[i]*v[i]; vPhi += phiq[i]*v[i];
} }
if (vPhi<=a)// on the left side of the node if (vPhi <= a)// on the left side of the node
{ {
if (y->nodeLeft() == nullptr)// left is a chemPoint if (y->nodeLeft() == nullptr)// left is a chemPoint
{ {
n2ndSearch_++; n2ndSearch_++;
if (y->leafLeft()->inEOA(phiq)) if (y->leafLeft()->inEOA(phiq))
{ {
x=y->leafLeft(); x = y->leafLeft();
return true; return true;
} }
} }
else// the left side is a node else // the left side is a node
{ {
if (inSubTree(phiq, y->nodeLeft(),x)) if (inSubTree(phiq, y->nodeLeft(),x))
{ {
@ -100,16 +100,16 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
// we reach the end of the subTree we can return the result // we reach the end of the subTree we can return the result
if (y->leafRight()->inEOA(phiq)) if (y->leafRight()->inEOA(phiq))
{ {
x=y->leafRight(); x = y->leafRight();
return true; return true;
} }
else else
{ {
x=nullptr; x = nullptr;
return false; return false;
} }
} }
else// test for n2ndSearch is done in the call of inSubTree else // test for n2ndSearch is done in the call of inSubTree
{ {
return inSubTree(phiq, y->nodeRight(),x); return inSubTree(phiq, y->nodeRight(),x);
} }
@ -124,11 +124,11 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
return true; return true;
} }
} }
else// the right side is a node else // the right side is a node
{ {
if (inSubTree(phiq, y->nodeRight(),x)) if (inSubTree(phiq, y->nodeRight(),x))
{ {
x=y->leafRight(); x = y->leafRight();
return true; return true;
} }
} }
@ -139,12 +139,12 @@ bool Foam::binaryTree<CompType, ThermoType>::inSubTree
n2ndSearch_++; n2ndSearch_++;
if (y->leafLeft()->inEOA(phiq)) if (y->leafLeft()->inEOA(phiq))
{ {
x=y->leafLeft(); x = y->leafLeft();
return true; return true;
} }
else else
{ {
x=nullptr; x = nullptr;
return false; return false;
} }
} }
@ -216,7 +216,7 @@ template<class CompType, class ThermoType>
Foam::chemPointISAT<CompType, ThermoType>* Foam::chemPointISAT<CompType, ThermoType>*
Foam::binaryTree<CompType, ThermoType>::chemPSibling(bn* y) Foam::binaryTree<CompType, ThermoType>::chemPSibling(bn* y)
{ {
if (y->parent()!=nullptr) if (y->parent() != nullptr)
{ {
if (y == y->parent()->nodeLeft())// y is on the left, return right side if (y == y->parent()->nodeLeft())// y is on the left, return right side
{ {
@ -235,6 +235,7 @@ Foam::binaryTree<CompType, ThermoType>::chemPSibling(bn* y)
return nullptr; return nullptr;
} }
} }
// the binaryNode is root_ and has no sibling // the binaryNode is root_ and has no sibling
return nullptr; return nullptr;
} }
@ -367,7 +368,9 @@ Foam::label Foam::binaryTree<CompType, ThermoType>::depth(bn* subTreeRoot)
} }
else else
{ {
return 1+max return
1
+ max
( (
depth(subTreeRoot->nodeLeft()), depth(subTreeRoot->nodeLeft()),
depth(subTreeRoot->nodeRight()) depth(subTreeRoot->nodeRight())
@ -379,14 +382,14 @@ Foam::label Foam::binaryTree<CompType, ThermoType>::depth(bn* subTreeRoot)
template<class CompType, class ThermoType> template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::insertNewLeaf void Foam::binaryTree<CompType, ThermoType>::insertNewLeaf
( (
const scalarField& phiq, const scalarField& phiq,
const scalarField& Rphiq, const scalarField& Rphiq,
const scalarSquareMatrix& A, const scalarSquareMatrix& A,
const scalarField& scaleFactor, const scalarField& scaleFactor,
const scalar& epsTol, const scalar& epsTol,
const label nCols, const label nCols,
chP*& phi0 chP*& phi0
) )
{ {
if (size_ == 0) // no points are stored if (size_ == 0) // no points are stored
{ {
@ -452,7 +455,7 @@ void Foam::binaryTree<CompType, ThermoType>::insertNewLeaf
root_ = newNode; root_ = newNode;
} }
phi0->node()=newNode; phi0->node() = newNode;
newChemPoint->node()=newNode; newChemPoint->node()=newNode;
} }
size_++; size_++;
@ -530,7 +533,7 @@ bool Foam::binaryTree<CompType, ThermoType>::secondaryBTSearch
n2ndSearch_++; n2ndSearch_++;
if (xS->inEOA(phiq)) if (xS->inEOA(phiq))
{ {
x=xS; x = xS;
return true; return true;
} }
} }
@ -557,7 +560,7 @@ bool Foam::binaryTree<CompType, ThermoType>::secondaryBTSearch
{ {
return true; return true;
} }
y=y->parent(); y = y->parent();
} }
// if we reach this point it is either because // if we reach this point it is either because
// we did not find another covering EOA in the entire tree or // we did not find another covering EOA in the entire tree or
@ -574,7 +577,6 @@ bool Foam::binaryTree<CompType, ThermoType>::secondaryBTSearch
template<class CompType, class ThermoType> template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::deleteLeaf(chP*& phi0) void Foam::binaryTree<CompType, ThermoType>::deleteLeaf(chP*& phi0)
{ {
if (size_ == 1) // only one point is stored if (size_ == 1) // only one point is stored
{ {
deleteDemandDrivenData(phi0); deleteDemandDrivenData(phi0);
@ -595,7 +597,7 @@ void Foam::binaryTree<CompType, ThermoType>::deleteLeaf(chP*& phi0)
root_->leafLeft()=siblingPhi0; root_->leafLeft()=siblingPhi0;
siblingPhi0->node()=root_; siblingPhi0->node()=root_;
} }
else if (z==z->parent()->nodeLeft()) else if (z == z->parent()->nodeLeft())
{ {
z->parent()->leafLeft() = siblingPhi0; z->parent()->leafLeft() = siblingPhi0;
z->parent()->nodeLeft() = nullptr; z->parent()->nodeLeft() = nullptr;
@ -642,16 +644,16 @@ void Foam::binaryTree<CompType, ThermoType>::balance()
//1) walk through the entire tree by starting with the tree's most left //1) walk through the entire tree by starting with the tree's most left
// chemPoint // chemPoint
chP* x=treeMin(); chP* x = treeMin();
List<chP*> chemPoints(size_); List<chP*> chemPoints(size_);
label chPi=0; label chPi=0;
//2) compute the mean composition //2) compute the mean composition
while(x!=nullptr) while (x != nullptr)
{ {
const scalarField& phij = x->phi(); const scalarField& phij = x->phi();
mean += phij; mean += phij;
chemPoints[chPi++] = x; chemPoints[chPi++] = x;
x=treeSuccessor(x); x = treeSuccessor(x);
} }
mean /= size_; mean /= size_;
@ -691,15 +693,15 @@ void Foam::binaryTree<CompType, ThermoType>::balance()
phiMaxDir.sort(); phiMaxDir.sort();
// delete reference to all node since the tree is reshaped // delete reference to all node since the tree is reshaped
deleteAllNode(); deleteAllNode();
root_=nullptr; root_ = nullptr;
// add the node for the two extremum // add the node for the two extremum
bn* newNode = new bn bn* newNode = new bn
( (
chemPoints[phiMaxDir.indices()[0]], chemPoints[phiMaxDir.indices()[0]],
chemPoints[phiMaxDir.indices()[phiMaxDir.size()-1]], chemPoints[phiMaxDir.indices()[phiMaxDir.size()-1]],
nullptr nullptr
); );
root_ = newNode; root_ = newNode;
chemPoints[phiMaxDir.indices()[0]]->node() = newNode; chemPoints[phiMaxDir.indices()[0]]->node() = newNode;
@ -719,8 +721,8 @@ void Foam::binaryTree<CompType, ThermoType>::balance()
new bn(phi0,chemPoints[phiMaxDir.indices()[cpi]], phi0->node()); new bn(phi0,chemPoints[phiMaxDir.indices()[cpi]], phi0->node());
// make the parent of phi0 point to the newly created node // make the parent of phi0 point to the newly created node
insertNode(phi0, nodeToAdd); insertNode(phi0, nodeToAdd);
phi0->node()=nodeToAdd; phi0->node() = nodeToAdd;
chemPoints[phiMaxDir.indices()[cpi]]->node()=nodeToAdd; chemPoints[phiMaxDir.indices()[cpi]]->node() = nodeToAdd;
} }
} }
@ -800,14 +802,14 @@ Foam::binaryTree<CompType, ThermoType>::treeSuccessor(chP* x)
template<class CompType, class ThermoType> template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::clear() void Foam::binaryTree<CompType, ThermoType>::clear()
{ {
// recursively delete the element in the subTree // Recursively delete the element in the subTree
deleteSubTree(); deleteSubTree();
// reset root node (should already be nullptr) // Reset root node (should already be nullptr)
root_=nullptr; root_ = nullptr;
// reset size_ // Reset size_
size_=0; size_ = 0;
} }
@ -818,4 +820,25 @@ bool Foam::binaryTree<CompType, ThermoType>::isFull()
} }
template<class CompType, class ThermoType>
void Foam::binaryTree<CompType, ThermoType>::resetNumRetrieve()
{
// Has to go along each chP of the tree
if (size_ > 0)
{
// First finds the first leaf
chP* chP0 = treeMin();
chP0->resetNumRetrieve();
// Then go to each successor
chP* nextChP = treeSuccessor(chP0);
while (nextChP != nullptr)
{
nextChP->resetNumRetrieve();
nextChP = treeSuccessor(nextChP);
}
}
}
// ************************************************************************* // // ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -230,6 +230,8 @@ public:
//- ListFull //- ListFull
bool isFull(); bool isFull();
void resetNumRetrieve();
}; };

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -33,7 +33,6 @@ License
template<class CompType, class ThermoType> template<class CompType, class ThermoType>
Foam::scalar Foam::chemPointISAT<CompType, ThermoType>::tolerance_; Foam::scalar Foam::chemPointISAT<CompType, ThermoType>::tolerance_;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class CompType, class ThermoType> template<class CompType, class ThermoType>
@ -75,7 +74,7 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
d[k] = -scale*sigma; d[k] = -scale*sigma;
for (label j=k+1; j<nCols; j++) for (label j=k+1; j<nCols; j++)
{ {
sum=0; sum = 0;
for ( label i=k; i<nCols; i++) for ( label i=k; i<nCols; i++)
{ {
sum += R(i, k)*R(i, j); sum += R(i, k)*R(i, j);
@ -137,11 +136,11 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrUpdate
} }
else if (mag(w[i]) > mag(w[i+1])) else if (mag(w[i]) > mag(w[i+1]))
{ {
w[i] = mag(w[i])*sqrt(1 + sqr(w[i+1]/w[i])); w[i] = mag(w[i])*sqrt(1.0 + sqr(w[i+1]/w[i]));
} }
else else
{ {
w[i] = mag(w[i+1])*sqrt(1 + sqr(w[i]/w[i+1])); w[i] = mag(w[i+1])*sqrt(1.0 + sqr(w[i]/w[i+1]));
} }
} }
@ -171,18 +170,18 @@ void Foam::chemPointISAT<CompType, ThermoType>::rotate
if (a == 0) if (a == 0)
{ {
c = 0; c = 0;
s = (b >= 0 ? 1.0 : -1.0); s = (b >= 0 ? 1 : -1);
} }
else if (mag(a) > mag(b)) else if (mag(a) > mag(b))
{ {
fact = b/a; fact = b/a;
c = sign(a)/sqrt(1 + sqr(fact)); c = sign(a)/sqrt(1.0 + sqr(fact));
s = fact*c; s = fact*c;
} }
else else
{ {
fact = a/b; fact = a/b;
s = sign(b)/sqrt(1 + sqr(fact)); s = sign(b)/sqrt(1.0 + sqr(fact));
c = fact*s; c = fact*s;
} }
for (label j=i;j<n;j++) for (label j=i;j<n;j++)
@ -220,20 +219,43 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
completeSpaceSize_(completeSpaceSize), completeSpaceSize_(completeSpaceSize),
nGrowth_(0), nGrowth_(0),
nActiveSpecies_(chemistry.mechRed()->NsSimp()), nActiveSpecies_(chemistry.mechRed()->NsSimp()),
completeToSimplifiedIndex_(completeSpaceSize-2),
simplifiedToCompleteIndex_(nActiveSpecies_), simplifiedToCompleteIndex_(nActiveSpecies_),
timeTag_(chemistry_.time().timeOutputValue()), timeTag_(chemistry_.timeSteps()),
lastTimeUsed_(chemistry_.time().timeOutputValue()), lastTimeUsed_(chemistry_.timeSteps()),
toRemove_(false), toRemove_(false),
maxNumNewDim_(coeffsDict.lookupOrDefault("maxNumNewDim",0)), maxNumNewDim_(coeffsDict.lookupOrDefault("maxNumNewDim",0)),
printProportion_(coeffsDict.lookupOrDefault("printProportion",false)) printProportion_(coeffsDict.lookupOrDefault("printProportion",false)),
numRetrieve_(0),
nLifeTime_(0),
variableTimeStep_
(
coeffsDict.lookupOrDefault("variableTimeStep", false)
),
completeToSimplifiedIndex_
(
completeSpaceSize - (2 + (variableTimeStep_ == 1 ? 1 : 0))
)
{ {
tolerance_=tolerance; tolerance_=tolerance;
if (this->variableTimeStep_)
{
nAdditionalEqns_ = 3;
iddeltaT_ = completeSpaceSize - 1;
scaleFactor_[iddeltaT_] *= phi_[iddeltaT_] / tolerance_;
}
else
{
nAdditionalEqns_ = 2;
iddeltaT_ = completeSpaceSize; // will not be used
}
idT_ = completeSpaceSize - nAdditionalEqns_;
idp_ = completeSpaceSize - nAdditionalEqns_ + 1;
bool isMechRedActive = chemistry_.mechRed()->active(); bool isMechRedActive = chemistry_.mechRed()->active();
if (isMechRedActive) if (isMechRedActive)
{ {
for (label i=0; i<completeSpaceSize-2; i++) for (label i=0; i<completeSpaceSize-nAdditionalEqns_; i++)
{ {
completeToSimplifiedIndex_[i] = completeToSimplifiedIndex_[i] =
chemistry.completeToSimplifiedIndex()[i]; chemistry.completeToSimplifiedIndex()[i];
@ -248,7 +270,7 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
label reduOrCompDim = completeSpaceSize; label reduOrCompDim = completeSpaceSize;
if (isMechRedActive) if (isMechRedActive)
{ {
reduOrCompDim = nActiveSpecies_ + 2; reduOrCompDim = nActiveSpecies_+nAdditionalEqns_;
} }
// SVD decomposition A = U*D*V^T // SVD decomposition A = U*D*V^T
@ -313,14 +335,32 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
completeSpaceSize_(p.completeSpaceSize()), completeSpaceSize_(p.completeSpaceSize()),
nGrowth_(p.nGrowth()), nGrowth_(p.nGrowth()),
nActiveSpecies_(p.nActiveSpecies()), nActiveSpecies_(p.nActiveSpecies()),
completeToSimplifiedIndex_(p.completeToSimplifiedIndex()),
simplifiedToCompleteIndex_(p.simplifiedToCompleteIndex()), simplifiedToCompleteIndex_(p.simplifiedToCompleteIndex()),
timeTag_(p.timeTag()), timeTag_(p.timeTag()),
lastTimeUsed_(p.lastTimeUsed()), lastTimeUsed_(p.lastTimeUsed()),
toRemove_(p.toRemove()), toRemove_(p.toRemove()),
maxNumNewDim_(p.maxNumNewDim()) maxNumNewDim_(p.maxNumNewDim()),
numRetrieve_(0),
nLifeTime_(0),
variableTimeStep_(p.variableTimeStep()),
completeToSimplifiedIndex_(p.completeToSimplifiedIndex())
{ {
tolerance_ = p.tolerance(); tolerance_ = p.tolerance();
if (this->variableTimeStep_)
{
nAdditionalEqns_ = 3;
idT_ = completeSpaceSize() - 3;
idp_ = completeSpaceSize() - 2;
iddeltaT_ = completeSpaceSize() - 1;
}
else
{
nAdditionalEqns_ = 2;
idT_ = completeSpaceSize() - 2;
idp_ = completeSpaceSize() - 1;
iddeltaT_ = completeSpaceSize(); // will not be used
}
} }
@ -331,11 +371,20 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
{ {
scalarField dphi(phiq-phi()); scalarField dphi(phiq-phi());
bool isMechRedActive = chemistry_.mechRed()->active(); bool isMechRedActive = chemistry_.mechRed()->active();
label dim = (isMechRedActive) ? nActiveSpecies_ : completeSpaceSize()-2; label dim(0);
scalar epsTemp=0; if (isMechRedActive)
{
dim = nActiveSpecies_;
}
else
{
dim = completeSpaceSize() - nAdditionalEqns_;
}
scalar epsTemp = 0;
List<scalar> propEps(completeSpaceSize(), scalar(0)); List<scalar> propEps(completeSpaceSize(), scalar(0));
for (label i=0; i<completeSpaceSize()-2; i++) for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
{ {
scalar temp = 0; scalar temp = 0;
@ -356,8 +405,12 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
temp += LT_(si, j)*dphi[sj]; temp += LT_(si, j)*dphi[sj];
} }
temp += LT_(si, nActiveSpecies_)*dphi[completeSpaceSize()-2]; temp += LT_(si, dim)*dphi[idT_];
temp += LT_(si, nActiveSpecies_+1)*dphi[completeSpaceSize()-1]; temp += LT_(si, dim+1)*dphi[idp_];
if (variableTimeStep_)
{
temp += LT_(si, dim+2)*dphi[iddeltaT_];
}
} }
else else
{ {
@ -373,28 +426,65 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
} }
// Temperature // Temperature
epsTemp += if (variableTimeStep_)
sqr {
( epsTemp +=
LT_(dim, dim)*dphi[completeSpaceSize()-2] sqr
+LT_(dim, dim+1)*dphi[completeSpaceSize()-1] (
); LT_(dim, dim)*dphi[idT_]
+LT_(dim, dim+1)*dphi[idp_]
+LT_(dim, dim+2)*dphi[iddeltaT_]
);
}
else
{
epsTemp +=
sqr
(
LT_(dim, dim)*dphi[idT_]
+LT_(dim, dim+1)*dphi[idp_]
);
}
// Pressure // Pressure
epsTemp += sqr(LT_(dim+1, dim+1)*dphi[completeSpaceSize()-1]); if (variableTimeStep_)
{
epsTemp +=
sqr
(
LT_(dim+1, dim+1)*dphi[idp_]
+LT_(dim+1, dim+2)*dphi[iddeltaT_]
);
}
else
{
epsTemp += sqr(LT_(dim+1, dim+1)*dphi[idp_]);
}
if (variableTimeStep_)
{
epsTemp += sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
}
if (printProportion_) if (printProportion_)
{ {
propEps[completeSpaceSize()-2] = propEps[idT_] = sqr
sqr
( (
LT_(dim, dim)*dphi[completeSpaceSize()-2] LT_(dim, dim)*dphi[idT_]
+ LT_(dim, dim+1)*dphi[completeSpaceSize()-1] + LT_(dim, dim+1)*dphi[idp_]
); );
propEps[completeSpaceSize()-1] = propEps[idp_] =
sqr(LT_(dim+1, dim+1)*dphi[completeSpaceSize()-1]); sqr(LT_(dim+1, dim+1)*dphi[idp_]);
if (variableTimeStep_)
{
propEps[iddeltaT_] =
sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
}
} }
if (sqrt(epsTemp) > 1 + tolerance_) if (sqrt(epsTemp) > 1 + tolerance_)
{ {
if (printProportion_) if (printProportion_)
@ -410,16 +500,20 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
} }
} }
word propName; word propName;
if (maxIndex >= completeSpaceSize()-2) if (maxIndex >= completeSpaceSize() - nAdditionalEqns_)
{ {
if(maxIndex == completeSpaceSize()-2) if (maxIndex == idT_)
{ {
propName = "T"; propName = "T";
} }
else if(maxIndex == completeSpaceSize()-1) else if (maxIndex == idp_)
{ {
propName = "p"; propName = "p";
} }
else if (maxIndex == iddeltaT_)
{
propName = "deltaT";
}
} }
else else
{ {
@ -461,7 +555,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
// Since we build only the solution for the species, T and p are not // Since we build only the solution for the species, T and p are not
// included // included
for (label i=0; i<completeSpaceSize()-2; i++) for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
{ {
dRl = 0; dRl = 0;
if (isMechRedActive) if (isMechRedActive)
@ -476,8 +570,12 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
label sj=simplifiedToCompleteIndex_[j]; label sj=simplifiedToCompleteIndex_[j];
dRl += Avar(si, j)*dphi[sj]; dRl += Avar(si, j)*dphi[sj];
} }
dRl += Avar(si, nActiveSpecies_)*dphi[completeSpaceSize()-2]; dRl += Avar(si, nActiveSpecies_)*dphi[idT_];
dRl += Avar(si, nActiveSpecies_+1)*dphi[completeSpaceSize()-1]; dRl += Avar(si, nActiveSpecies_+1)*dphi[idp_];
if (variableTimeStep_)
{
dRl += Avar(si, nActiveSpecies_+2)*dphi[iddeltaT_];
}
} }
else else
{ {
@ -522,7 +620,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
// check if the difference of active species is lower than the maximum // check if the difference of active species is lower than the maximum
// number of new dimensions allowed // number of new dimensions allowed
for (label i=0; i<completeSpaceSize()-2; i++) for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
{ {
// first test if the current chemPoint has an inactive species // first test if the current chemPoint has an inactive species
// corresponding to an active one in the query point // corresponding to an active one in the query point
@ -588,8 +686,8 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
{ {
scalarSquareMatrix LTvar = LT_; // take a copy of LT_ scalarSquareMatrix LTvar = LT_; // take a copy of LT_
scalarSquareMatrix Avar = A_; // take a copy of A_ scalarSquareMatrix Avar = A_; // take a copy of A_
LT_ = scalarSquareMatrix(nActiveSpecies_+2, Zero); LT_ = scalarSquareMatrix(nActiveSpecies_+nAdditionalEqns_, Zero);
A_ = scalarSquareMatrix(nActiveSpecies_+2, Zero); A_ = scalarSquareMatrix(nActiveSpecies_+nAdditionalEqns_, Zero);
// write the initial active species // write the initial active species
for (label i=0; i<initNActiveSpecies; i++) for (label i=0; i<initNActiveSpecies; i++)
@ -621,6 +719,13 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
LTvar(initNActiveSpecies+1, initNActiveSpecies+1); LTvar(initNActiveSpecies+1, initNActiveSpecies+1);
A_(nActiveSpecies_+1, nActiveSpecies_+1)= A_(nActiveSpecies_+1, nActiveSpecies_+1)=
Avar(initNActiveSpecies+1, initNActiveSpecies+1); Avar(initNActiveSpecies+1, initNActiveSpecies+1);
if (variableTimeStep_)
{
LT_(nActiveSpecies_+2, nActiveSpecies_+2)=
LTvar(initNActiveSpecies+2, initNActiveSpecies+2);
A_(nActiveSpecies_+2, nActiveSpecies_+2)=
Avar(initNActiveSpecies+2, initNActiveSpecies+2);
}
for (label i=initNActiveSpecies; i<nActiveSpecies_;i++) for (label i=initNActiveSpecies; i<nActiveSpecies_;i++)
{ {
@ -631,7 +736,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
} }
} }
dim = nActiveSpecies_ + 2; dim = nActiveSpecies_ + nAdditionalEqns_;
} }
// beginning of grow algorithm // beginning of grow algorithm
@ -641,7 +746,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
for (label i=0; i<dim; i++) for (label i=0; i<dim; i++)
{ {
for (label j=i; j<dim-2; j++)// LT is upper triangular for (label j=i; j<dim-nAdditionalEqns_; j++)// LT is upper triangular
{ {
label sj = j; label sj = j;
if (isMechRedActive) if (isMechRedActive)
@ -650,8 +755,12 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
} }
phiTilde[i] += LT_(i, j)*dphi[sj]; phiTilde[i] += LT_(i, j)*dphi[sj];
} }
phiTilde[i] += LT_(i, dim-2)*dphi[completeSpaceSize()-2]; phiTilde[i] += LT_(i, dim-nAdditionalEqns_)*dphi[idT_];
phiTilde[i] += LT_(i, dim-1)*dphi[completeSpaceSize()-1]; phiTilde[i] += LT_(i, dim-nAdditionalEqns_+1)*dphi[idp_];
if (variableTimeStep_)
{
phiTilde[i] += LT_(i, dim-nAdditionalEqns_ + 2)*dphi[iddeltaT_];
}
normPhiTilde += sqr(phiTilde[i]); normPhiTilde += sqr(phiTilde[i]);
} }
@ -678,4 +787,55 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
} }
template<class CompType, class ThermoType>
void Foam::chemPointISAT<CompType, ThermoType>::increaseNumRetrieve()
{
this->numRetrieve_++;
}
template<class CompType, class ThermoType>
void Foam::chemPointISAT<CompType, ThermoType>::resetNumRetrieve()
{
this->numRetrieve_ = 0;
}
template<class CompType, class ThermoType>
void Foam::chemPointISAT<CompType, ThermoType>::increaseNLifeTime()
{
this->nLifeTime_++;
}
template<class CompType, class ThermoType>
Foam::label Foam::chemPointISAT<CompType, ThermoType>::
simplifiedToCompleteIndex
(
const label i
)
{
if (i < nActiveSpecies_)
{
return simplifiedToCompleteIndex_[i];
}
else if (i == nActiveSpecies_)
{
return completeSpaceSize_-nAdditionalEqns_;
}
else if (i == nActiveSpecies_ + 1)
{
return completeSpaceSize_-nAdditionalEqns_ + 1;
}
else if (variableTimeStep_ && (i == nActiveSpecies_ + 2))
{
return completeSpaceSize_-nAdditionalEqns_ + 2;
}
else
{
return -1;
}
}
// ************************************************************************* // // ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -146,6 +146,7 @@ class chemPointISAT
TDACChemistryModel<CompType, ThermoType>& chemistry_; TDACChemistryModel<CompType, ThermoType>& chemistry_;
//- Vector storing the composition, temperature and pressure //- Vector storing the composition, temperature and pressure
// and deltaT if a variable time step is set on
scalarField phi_; scalarField phi_;
//- Vector storing the mapping of the composition phi //- Vector storing the mapping of the composition phi
@ -176,14 +177,12 @@ class chemPointISAT
//- Number of active species stored in the chemPoint //- Number of active species stored in the chemPoint
label nActiveSpecies_; label nActiveSpecies_;
//- Vectors that store the index conversion between the simplified
// Vectors that store the index conversion between the simplified // and the complete chemical mechanism
// and the complete chemical mechanism
List<label> completeToSimplifiedIndex_;
List<label> simplifiedToCompleteIndex_; List<label> simplifiedToCompleteIndex_;
scalar timeTag_; label timeTag_;
scalar lastTimeUsed_; label lastTimeUsed_;
bool toRemove_; bool toRemove_;
@ -191,6 +190,26 @@ class chemPointISAT
Switch printProportion_; Switch printProportion_;
//- Variable to store the number of retrieves the chemPoint
// will generate at each time step
label numRetrieve_;
//- Variable to store the number of time steps the chempoint is allowed
// to still live according to the maxChPLifeTime_ parameter
label nLifeTime_;
//- Switch to allow variable time step (off by default)
Switch variableTimeStep_;
List<label> completeToSimplifiedIndex_;
//- Number of equations in addition to the species eqs.
label nAdditionalEqns_;
label idT_;
label idp_;
label iddeltaT_;
//- QR decomposition of a matrix //- QR decomposition of a matrix
// Input : nCols cols number // Input : nCols cols number
// R the matrix to decompose // R the matrix to decompose
@ -311,6 +330,7 @@ public:
{ {
return A_; return A_;
} }
inline const scalarSquareMatrix& LT() const inline const scalarSquareMatrix& LT() const
{ {
return LT_; return LT_;
@ -336,32 +356,23 @@ public:
return simplifiedToCompleteIndex_; return simplifiedToCompleteIndex_;
} }
inline label simplifiedToCompleteIndex(label i) //- Increases the number of retrieves the chempoint has generated
{ void increaseNumRetrieve();
if (i < nActiveSpecies_)
{
return simplifiedToCompleteIndex_[i];
}
else if (i == nActiveSpecies_)
{
return completeSpaceSize_-2;
}
else if (i == nActiveSpecies_+1)
{
return completeSpaceSize_-1;
}
else
{
return -1;
}
}
inline const scalar& timeTag() //- Resets the number of retrieves at each time step
void resetNumRetrieve();
//- Increases the "counter" of the chP life
void increaseNLifeTime();
label simplifiedToCompleteIndex(const label i);
inline const label& timeTag()
{ {
return timeTag_; return timeTag_;
} }
inline scalar& lastTimeUsed() inline label& lastTimeUsed()
{ {
return lastTimeUsed_; return lastTimeUsed_;
} }
@ -376,6 +387,21 @@ public:
return maxNumNewDim_; return maxNumNewDim_;
} }
inline const label& numRetrieve()
{
return numRetrieve_;
}
inline const label& nLifeTime()
{
return nLifeTime_;
}
inline Switch variableTimeStep()
{
return variableTimeStep_;
}
// ISAT functions // ISAT functions
//- To RETRIEVE the mapping from the stored chemPoint phi, the query //- To RETRIEVE the mapping from the stored chemPoint phi, the query

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -40,7 +40,8 @@ Foam::chemistryTabulationMethod<CompType, ThermoType>::chemistryTabulationMethod
active_(coeffsDict_.lookupOrDefault<Switch>("active", false)), active_(coeffsDict_.lookupOrDefault<Switch>("active", false)),
log_(coeffsDict_.lookupOrDefault<Switch>("log", false)), log_(coeffsDict_.lookupOrDefault<Switch>("log", false)),
chemistry_(chemistry), chemistry_(chemistry),
tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4)) tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4)),
variableTimeStep_(coeffsDict_.lookupOrDefault("variableTimeStep", false))
{} {}

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -74,6 +74,8 @@ protected:
scalar tolerance_; scalar tolerance_;
//- Switch to allow variable time step (off by default)
const Switch variableTimeStep_;
public: public:
@ -130,6 +132,11 @@ public:
return active_ && log_; return active_ && log_;
} }
inline bool variableTimeStep()
{
return variableTimeStep_;
}
inline scalar tolerance() const inline scalar tolerance() const
{ {
return tolerance_; return tolerance_;
@ -152,11 +159,12 @@ public:
// Add function: (only virtual here) // Add function: (only virtual here)
// Add information to the tabulation algorithm. Give the reference for // Add information to the tabulation algorithm. Give the reference for
// future retrieve (phiQ) and the corresponding result (RphiQ). // future retrieve (phiQ) and the corresponding result (RphiQ).
virtual bool add virtual label add
( (
const scalarField& phiQ, const scalarField& phiQ,
const scalarField& RphiQ, const scalarField& RphiQ,
const scalar rho const scalar rho,
const scalar deltaT
) = 0; ) = 0;
// Update function: (only virtual here) // Update function: (only virtual here)

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -104,11 +104,12 @@ public:
// existing point or add a new leaf to the binary tree Input : phiq // existing point or add a new leaf to the binary tree Input : phiq
// the new composition to store Rphiq the mapping of the new // the new composition to store Rphiq the mapping of the new
// composition point // composition point
virtual bool add virtual label add
( (
const scalarField& phiq, const scalarField& phiq,
const scalarField& Rphiq, const scalarField& Rphiq,
const scalar rho const scalar rho,
const scalar deltaT
) )
{ {
NotImplemented; NotImplemented;

View File

@ -0,0 +1,48 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object CH4;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0.0;
boundaryField
{
fuel
{
type fixedValue;
value uniform 1.0;
}
air
{
type fixedValue;
value uniform 0.0;
}
outlet
{
type inletOutlet;
inletValue uniform 0.0;
value uniform 0.0;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,48 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object CO2;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
fuel
{
type fixedValue;
value uniform 0;
}
air
{
type fixedValue;
value uniform 0;
}
outlet
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,48 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object H2O;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
fuel
{
type fixedValue;
value uniform 0;
}
air
{
type fixedValue;
value uniform 0;
}
outlet
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,48 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object N2;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 1;
boundaryField
{
fuel
{
type fixedValue;
value uniform 0;
}
air
{
type fixedValue;
value uniform 0.77;
}
outlet
{
type inletOutlet;
inletValue uniform 1;
value uniform 1;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,47 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object O2;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
fuel
{
type fixedValue;
value uniform 0.0;
}
air
{
type fixedValue;
value uniform 0.23;
}
outlet
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,47 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0 0 0];
internalField uniform 2000;
boundaryField
{
fuel
{
type fixedValue;
value uniform 293;
}
air
{
type fixedValue;
value uniform 293;
}
outlet
{
type inletOutlet;
inletValue uniform 293;
value uniform 293;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,46 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
fuel
{
type fixedValue;
value uniform (0.1 0 0);
}
air
{
type fixedValue;
value uniform (-0.1 0 0);
}
outlet
{
type pressureInletOutletVelocity;
value $internalField;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,48 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object Ydefault;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0.0;
boundaryField
{
fuel
{
type fixedValue;
value uniform 0.0;
}
air
{
type fixedValue;
value uniform 0.0;
}
outlet
{
type inletOutlet;
inletValue uniform 0.0;
value uniform 0.0;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,45 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object alphat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
fuel
{
type fixedValue;
value uniform 0;
}
air
{
type fixedValue;
value uniform 0;
}
outlet
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,44 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 1e5;
boundaryField
{
fuel
{
type zeroGradient;
}
air
{
type zeroGradient;
}
outlet
{
type totalPressure;
p0 $internalField;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,137 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object chemistryProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
chemistryType
{
chemistrySolver ode;
chemistryThermo psi;
TDAC on;
}
chemistry on;
initialChemicalTimeStep 1e-7;
odeCoeffs
{
solver seulex;
absTol 1e-12;
relTol 1e-1;
}
reduction
{
// Activate reduction
active on;
// Switch logging of the reduction statistics and performance
log on;
// Tolerance depends on the reduction method, see details for each method
tolerance 1e-4;
// Available methods: DRG, DAC, DRGEP, PFA, EFA
method DAC;
// Search initiating set (SIS) of species, needed for most methods
initialSet
{
CO;
CH4;
HO2;
}
// For DAC, option to automatically change the SIS switch from HO2 to H2O
// and CO to CO2, + disable fuel
automaticSIS off;
// When automaticSIS, the method needs to know the fuel
fuelSpecies
{
CH4 1;
}
}
tabulation
{
// Activate tabulation
active on;
// Switch logging of the tabulation statistics and performance
log on;
printProportion off;
printNumRetrieve off;
variableTimeStep on;
// Tolerance used for retrieve and grow
tolerance 1e-3;
// ISAT is the only method currently available
method ISAT;
// Scale factors used in the definition of the ellipsoid of accuracy
scaleFactor
{
otherSpecies 1;
Temperature 25000;
Pressure 1e15;
deltaT 1;
}
// Maximum number of leafs stored in the binary tree
maxNLeafs 2000;
// Maximum life time of the leafs (in time steps) used in unsteady
// simulations to force renewal of the stored chemPoints and keep the tree
// small
chPMaxLifeTime 100;
// Maximum number of growth allowed on a chemPoint to avoid distorted
// chemPoints
maxGrowth 10;
// Number of time steps between analysis of the tree to remove old
// chemPoints or try to balance it
checkEntireTreeInterval 5;
// Parameters used to decide whether to balance or not if the tree's depth
// is larger than maxDepthFactor*log2(nLeafs) then balance the tree
maxDepthFactor 2;
// Try to balance the tree only if the size of the tree is greater
minBalanceThreshold 30;
// Activate the use of a MRU (most recently used) list
MRURetrieve false;
// Maximum size of the MRU list
maxMRUSize 0;
// Allow to grow points
growPoints true;
// When mechanism reduction is used, new dimensions might be added
// maxNumNewDim set the maximum number of new dimensions added during a
// growth
maxNumNewDim 10;
}
// ************************************************************************* //

View File

@ -0,0 +1,27 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object combustionProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
combustionModel laminar<psiChemistryCombustion>;
active true;
laminarCoeffs
{
}
// ************************************************************************* //

View File

@ -0,0 +1,28 @@
elements
(
O
C
H
N
);
species
(
O2
H2O
CH4
CO2
N2
);
reactions
{
methaneReaction
{
type irreversibleArrheniusReaction;
reaction "CH4 + 2O2 = CO2 + 2H2O";
A 5.2e16;
beta 0;
Ta 14906;
}
}

View File

@ -0,0 +1,152 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 3.0.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object thermo.compressibleGas;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
O2
{
specie
{
nMoles 1;
molWeight 31.9988;
}
elements
{
O 2;
}
thermodynamics
{
Tlow 200;
Thigh 5000;
Tcommon 1000;
highCpCoeffs ( 3.69758 0.00061352 -1.25884e-07 1.77528e-11 -1.13644e-15 -1233.93 3.18917 );
lowCpCoeffs ( 3.21294 0.00112749 -5.75615e-07 1.31388e-09 -8.76855e-13 -1005.25 6.03474 );
}
transport
{
As 1.67212e-06;
Ts 170.672;
}
}
H2O
{
specie
{
nMoles 1;
molWeight 18.0153;
}
elements
{
O 1;
H 2;
}
thermodynamics
{
Tlow 200;
Thigh 5000;
Tcommon 1000;
highCpCoeffs ( 2.67215 0.00305629 -8.73026e-07 1.201e-10 -6.39162e-15 -29899.2 6.86282 );
lowCpCoeffs ( 3.38684 0.00347498 -6.3547e-06 6.96858e-09 -2.50659e-12 -30208.1 2.59023 );
}
transport
{
As 1.67212e-06;
Ts 170.672;
}
}
CH4
{
specie
{
nMoles 1;
molWeight 16.0428;
}
elements
{
C 1;
H 4;
}
thermodynamics
{
Tlow 200;
Thigh 6000;
Tcommon 1000;
highCpCoeffs ( 1.63543 0.0100844 -3.36924e-06 5.34973e-10 -3.15528e-14 -10005.6 9.9937 );
lowCpCoeffs ( 5.14988 -0.013671 4.91801e-05 -4.84744e-08 1.66694e-11 -10246.6 -4.64132 );
}
transport
{
As 1.67212e-06;
Ts 170.672;
}
}
CO2
{
specie
{
nMoles 1;
molWeight 44.01;
}
elements
{
C 1;
O 2;
}
thermodynamics
{
Tlow 200;
Thigh 5000;
Tcommon 1000;
highCpCoeffs ( 4.45362 0.00314017 -1.27841e-06 2.394e-10 -1.66903e-14 -48967 -0.955396 );
lowCpCoeffs ( 2.27572 0.00992207 -1.04091e-05 6.86669e-09 -2.11728e-12 -48373.1 10.1885 );
}
transport
{
As 1.67212e-06;
Ts 170.672;
}
}
N2
{
specie
{
nMoles 1;
molWeight 28.0134;
}
elements
{
N 2;
}
thermodynamics
{
Tlow 200;
Thigh 5000;
Tcommon 1000;
highCpCoeffs ( 2.92664 0.00148798 -5.68476e-07 1.0097e-10 -6.75335e-15 -922.798 5.98053 );
lowCpCoeffs ( 3.29868 0.00140824 -3.96322e-06 5.64152e-09 -2.44486e-12 -1020.9 3.95037 );
}
transport
{
As 1.67212e-06;
Ts 170.672;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,36 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType
{
type hePsiThermo;
mixture reactingMixture;
transport sutherland;
thermo janaf;
energy sensibleEnthalpy;
equationOfState perfectGas;
specie specie;
}
inertSpecie N2;
chemistryReader foamChemistryReader;
foamChemistryFile "$FOAM_CASE/constant/reactionsGRI";
foamChemistryThermoFile "$FOAM_CASE/constant/thermo.compressibleGasGRI";
// ************************************************************************* //

View File

@ -0,0 +1,21 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType laminar;
// ************************************************************************* //

View File

@ -0,0 +1,82 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.0 -0.01 -0.01)
(0.02 -0.01 -0.01)
(0.02 0.01 -0.01)
(0.0 0.01 -0.01)
(0.0 -0.01 0.01)
(0.02 -0.01 0.01)
(0.02 0.01 0.01)
(0.0 0.01 0.01)
);
blocks
(
hex (0 1 2 3 4 5 6 7) (100 40 1) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
fuel
{
type patch;
faces
(
(0 4 7 3)
);
}
air
{
type patch;
faces
(
(1 2 6 5)
);
}
outlet
{
type patch;
faces
(
(0 1 5 4)
(7 6 2 3)
);
}
frontAndBack
{
type empty;
faces
(
(4 5 6 7)
(0 3 2 1)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //

View File

@ -0,0 +1,48 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application reactingFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 1000;
deltaT 1;
writeControl runTime;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //

View File

@ -0,0 +1,29 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
note "mesh decomposition control dictionary";
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 4;
method hierarchical;
hierarchicalCoeffs
{
n (2 2 1);
delta 0.001;
order xyz;
}
// ************************************************************************* //

View File

@ -0,0 +1,57 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default localEuler;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss limitedLinearV 1;
div(phi,Yi_h) Gauss limitedLinear 1;
div(phi,K) Gauss limitedLinear 1;
div(phid,p) Gauss limitedLinear 1;
div(phi,epsilon) Gauss limitedLinear 1;
div(phi,k) Gauss limitedLinear 1;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //

View File

@ -0,0 +1,83 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
"rho.*"
{
solver diagonal;
}
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-6;
relTol 0.1;
}
pFinal
{
$p;
relTol 0;
}
"(U|h|k|epsilon)"
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-6;
relTol 0.1;
}
"(U|h|k|epsilon)Final"
{
$U;
relTol 0.1;
}
Yi
{
$U;
relTol 0.1;
}
}
PIMPLE
{
momentumPredictor no;
nOuterCorrectors 1;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
maxDeltaT 1e-4;
maxCo 1;
alphaTemp 0.05;
rDeltaTSmoothingCoeff 1;
rDeltaTDampingCoeff 1;
}
relaxationFactors
{
equations
{
".*" 1;
}
}
// ************************************************************************* //