zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

chemistryModel: Refactored interface to chemistry reduction

Browse Source 
in preparation for un-templating ISAT
This commit is contained in:
Henry Weller
2022-01-18 22:06:04 +00:00
parent 427307c190
commit a3f7c38835
12 changed files with 107 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
View File

@ -867,7 +867,14 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
Rphiq[Rphiq.size()-2] = p;
Rphiq[Rphiq.size()-1] = deltaT[celli];
tabulation_.add(phiq, Rphiq, celli, rho0, deltaT[celli]);
tabulation_.add
(
phiq,
Rphiq,
mechRed_.nActiveSpecies(),
celli,
deltaT[celli]
);
}
// When operations are done and if mechanism reduction is active,

5
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H
View File

@ -330,9 +330,8 @@ public:
// Mechanism reduction access functions
//- Return access to the mechanism reduction method
// Needed by ISAT
inline const chemistryReductionMethod<ThermoType>& mechRed() const;
//- Return true if mechanism reduction is active
inline bool reduction() const;
//- Allow the reduction method to reset the number of species
inline void setNSpecie(const label newNs);

7
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModelI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -119,10 +119,9 @@ Foam::chemistryModel<ThermoType>::Y() const
template<class ThermoType>
inline const Foam::chemistryReductionMethod<ThermoType>&
Foam::chemistryModel<ThermoType>::mechRed() const
inline bool Foam::chemistryModel<ThermoType>::reduction() const
{
return mechRed_;
return mechRedActive_;
}

11
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.C
View File

@ -43,6 +43,7 @@ Foam::chemistryTabulationMethods::ISAT<ThermoType>::ISAT
coeffsDict_(chemistryProperties.subDict("tabulation")),
chemistry_(chemistry),
log_(coeffsDict_.lookupOrDefault<Switch>("log", false)),
reduction_(chemistry_.reduction()),
chemisTree_(*this, coeffsDict_),
scaleFactor_(chemistry.nEqns() + 1, 1),
runTime_(chemistry.time()),
@ -198,7 +199,6 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::calcNewC
)
{
const label nEqns = chemistry_.nEqns(); // Species, T, p
const bool mechRedActive = chemistry_.mechRed().active();
const List<label>& completeToSimplified = phi0->completeToSimplifiedIndex();
const scalarField dphi(phiq - phi0->phi());
@ -213,12 +213,12 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::calcNewC
Rphiq = phi0->Rphi();
for (label i=0; i<nEqns + 1; i++)
{
if (mechRedActive)
if (reduction_)
{
const label si =
i < nEqns - 2
? completeToSimplified[i]
: i - (nEqns - 2) + phi0->nActiveSpecies();
: i - (nEqns - 2) + phi0->nActive();
if (si != -1)
{
@ -229,7 +229,7 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::calcNewC
const label sj =
j < nEqns - 2
? completeToSimplified[j]
: j - (nEqns - 2) + phi0->nActiveSpecies();
: j - (nEqns - 2) + phi0->nActive();
if (sj != -1)
{
@ -490,6 +490,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
(
const scalarField& phiq,
const scalarField& Rphiq,
const label nActive,
const label li,
const scalar deltaT
)
@ -568,6 +569,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
scaleFactor(),
tolerance_,
scaleFactor_.size(),
nActive,
nulPhi
);
deleteDemandDrivenData(tempList[i]);
@ -592,6 +594,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
scaleFactor(),
tolerance_,
scaleFactor_.size(),
nActive,
lastSearch_ // lastSearch_ may be nullptr (handled by binaryTree)
);
if (lastSearch_ != nullptr)

10
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.H
View File

@ -68,6 +68,9 @@ class ISAT
//- Switch to select performance logging
Switch log_;
//- Is reduction applied to the state vectors
const bool reduction_;
//- List of the stored 'points' organised in a binary tree
binaryTree<ThermoType> chemisTree_;
@ -232,6 +235,12 @@ public:
return true;
}
//- Return true if reduction is applied to the state variables
bool reduction() const
{
return reduction_;
}
const chemistryModel<ThermoType>& chemistry()
{
return chemistry_;
@ -271,6 +280,7 @@ public:
(
const scalarField& phiq,
const scalarField& Rphiq,
const label nActive,
const label li,
const scalar deltaT
);

22
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/binaryNode/binaryNode.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -70,7 +70,7 @@ void Foam::binaryNode<ThermoType>::calcV
{
// LT is the transpose of the L matrix
scalarSquareMatrix& LT = elementLeft->LT();
bool mechReductionActive = elementLeft->chemistry().mechRed().active();
const bool reduction = elementLeft->table().reduction();
// Difference of composition in the full species domain
scalarField phiDif(elementRight->phi() - elementLeft->phi());
@ -82,7 +82,7 @@ void Foam::binaryNode<ThermoType>::calcV
{
label si = i;
bool outOfIndexI = true;
if (mechReductionActive)
if (reduction)
{
if (i<elementLeft->completeSpaceSize() - 3)
{
@ -93,17 +93,18 @@ void Foam::binaryNode<ThermoType>::calcV
{
outOfIndexI = false;
const label dif = i - (elementLeft->completeSpaceSize() - 3);
si = elementLeft->nActiveSpecies() + dif;
si = elementLeft->nActive() + dif;
}
}
if (!mechReductionActive || (mechReductionActive && !(outOfIndexI)))
if (!reduction || (reduction && !(outOfIndexI)))
{
v[i] = 0;
for (label j=0; j<elementLeft->completeSpaceSize(); j++)
{
label sj = j;
bool outOfIndexJ = true;
if (mechReductionActive)
if (reduction)
{
if (j < elementLeft->completeSpaceSize() - 3)
{
@ -115,14 +116,11 @@ void Foam::binaryNode<ThermoType>::calcV
outOfIndexJ = false;
const label dif =
j - (elementLeft->completeSpaceSize() - 3);
sj = elementLeft->nActiveSpecies() + dif;
sj = elementLeft->nActive() + dif;
}
}
if
(
!mechReductionActive
||(mechReductionActive && !(outOfIndexJ))
)
if (!reduction ||(reduction && !(outOfIndexJ)))
{
// Since L is a lower triangular matrix k=0->min(i, j)
for (label k=0; k<=min(si, sj); k++)

5
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/binaryTree/binaryTree.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -382,6 +382,7 @@ void Foam::binaryTree<ThermoType>::insertNewLeaf
const scalarField& scaleFactor,
const scalar& epsTol,
const label nCols,
const label nActive,
chP*& phi0
)
{
@ -401,6 +402,7 @@ void Foam::binaryTree<ThermoType>::insertNewLeaf
scaleFactor,
epsTol,
nCols,
nActive,
coeffsDict_,
root_
);
@ -428,6 +430,7 @@ void Foam::binaryTree<ThermoType>::insertNewLeaf
scaleFactor,
epsTol,
nCols,
nActive,
coeffsDict_
);
// insert new node on the parent node in the position of the

3
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/binaryTree/binaryTree.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -188,6 +188,7 @@ public:
const scalarField& scaleFactor,
const scalar& epsTol,
const label nCols,
const label nActive,
chP*& phi0
);

95
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -204,8 +204,9 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
const scalarField& Rphi,
const scalarSquareMatrix& A,
const scalarField& scaleFactor,
const scalar& tolerance,
const label& completeSpaceSize,
const scalar tolerance,
const label completeSpaceSize,
const label nActive,
const dictionary& coeffsDict,
binaryNode<ThermoType>* node
)
@ -218,8 +219,8 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
node_(node),
completeSpaceSize_(completeSpaceSize),
nGrowth_(0),
nActiveSpecies_(table_.chemistry().mechRed().nActiveSpecies()),
simplifiedToCompleteIndex_(nActiveSpecies_),
nActive_(nActive),
simplifiedToCompleteIndex_(nActive_),
timeTag_(table.timeSteps()),
lastTimeUsed_(table.timeSteps()),
toRemove_(false),
@ -237,21 +238,20 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
idT_ = completeSpaceSize - 3;
idp_ = completeSpaceSize - 2;
const bool mechRedActive = table_.chemistry().mechRed().active();
if (mechRedActive)
if (table_.reduction())
{
for (label i=0; i<completeSpaceSize-3; i++)
{
completeToSimplifiedIndex_[i] = table_.chemistry().cTos(i);
}
for (label i=0; i<nActiveSpecies_; i++)
for (label i=0; i<nActive_; i++)
{
simplifiedToCompleteIndex_[i] = table_.chemistry().sToc(i);
}
}
const label reduOrCompDim =
mechRedActive ? nActiveSpecies_ + 3 : completeSpaceSize;
table_.reduction() ? nActive_ + 3 : completeSpaceSize;
// SVD decomposition A = U*D*V^T
SVD svdA(A);
@ -277,7 +277,7 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
{
label compi = i;
if (mechRedActive)
if (table_.reduction())
{
compi = simplifiedToCompleteIndex(i);
}
@ -314,7 +314,7 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
node_(p.node()),
completeSpaceSize_(p.completeSpaceSize()),
nGrowth_(p.nGrowth()),
nActiveSpecies_(p.nActiveSpecies()),
nActive_(p.nActive()),
simplifiedToCompleteIndex_(p.simplifiedToCompleteIndex()),
timeTag_(p.timeTag()),
lastTimeUsed_(p.lastTimeUsed()),
@ -338,10 +338,9 @@ template<class ThermoType>
bool Foam::chemPointISAT<ThermoType>::inEOA(const scalarField& phiq)
{
const scalarField dphi(phiq - phi());
const bool mechRedActive = table_.chemistry().mechRed().active();
const label dim =
mechRedActive ? nActiveSpecies_ : completeSpaceSize() - 3;
table_.reduction() ? nActive_ : completeSpaceSize() - 3;
scalar epsTemp = 0;
List<scalar> propEps(completeSpaceSize(), scalar(0));
@ -355,19 +354,19 @@ bool Foam::chemPointISAT<ThermoType>::inEOA(const scalarField& phiq)
// inactive species), just multiply the diagonal element and dphi
if
(
!(mechRedActive)
||(mechRedActive && completeToSimplifiedIndex_[i] != -1)
!(table_.reduction())
||(table_.reduction() && completeToSimplifiedIndex_[i] != -1)
)
{
const label si =
mechRedActive
table_.reduction()
? completeToSimplifiedIndex_[i]
: i;
for (label j=si; j<dim; j++)// LT is upper triangular
{
const label sj =
mechRedActive
table_.reduction()
? simplifiedToCompleteIndex_[j]
: j;
@ -487,18 +486,17 @@ bool Foam::chemPointISAT<ThermoType>::checkSolution
const scalarField dphi(phiq - phi());
const scalarField& scaleFactorV(scaleFactor());
const scalarSquareMatrix& Avar(A());
const bool mechRedActive = table_.chemistry().mechRed().active();
scalar dRl = 0;
const label dim =
mechRedActive ? nActiveSpecies_ : completeSpaceSize() - 2;
table_.reduction() ? nActive_ : completeSpaceSize() - 2;
// Since we build only the solution for the species, T and p are not
// included
for (label i=0; i<completeSpaceSize()-3; i++)
{
dRl = 0;
if (mechRedActive)
if (table_.reduction())
{
const label si = completeToSimplifiedIndex_[i];
@ -510,9 +508,9 @@ bool Foam::chemPointISAT<ThermoType>::checkSolution
const label sj = simplifiedToCompleteIndex_[j];
dRl += Avar(si, j)*dphi[sj];
}
dRl += Avar(si, nActiveSpecies_)*dphi[idT_];
dRl += Avar(si, nActiveSpecies_+1)*dphi[idp_];
dRl += Avar(si, nActiveSpecies_+2)*dphi[iddeltaT_];
dRl += Avar(si, nActive_)*dphi[idT_];
dRl += Avar(si, nActive_+1)*dphi[idp_];
dRl += Avar(si, nActive_+2)*dphi[iddeltaT_];
}
else
{
@ -546,10 +544,9 @@ template<class ThermoType>
bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
{
const scalarField dphi(phiq - phi());
const label initNActiveSpecies(nActiveSpecies_);
const bool mechRedActive = table_.chemistry().mechRed().active();
const label initNActiveSpecies(nActive_);
if (mechRedActive)
if (table_.reduction())
{
label activeAdded(0);
DynamicList<label> dimToAdd(0);
@ -602,11 +599,11 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
}
// the number of added dimension to the current chemPoint
nActiveSpecies_ += dimToAdd.size();
simplifiedToCompleteIndex_.setSize(nActiveSpecies_);
nActive_ += dimToAdd.size();
simplifiedToCompleteIndex_.setSize(nActive_);
forAll(dimToAdd, i)
{
label si = nActiveSpecies_ - dimToAdd.size() + i;
label si = nActive_ - dimToAdd.size() + i;
// add the new active species
simplifiedToCompleteIndex_[si] = dimToAdd[i];
completeToSimplifiedIndex_[dimToAdd[i]] = si;
@ -618,12 +615,12 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
// (change the diagonal position)
//-set all element of the new lines and columns to zero except diagonal
// (=1/(tolerance*scaleFactor))
if (nActiveSpecies_ > initNActiveSpecies)
if (nActive_ > initNActiveSpecies)
{
const scalarSquareMatrix LTvar = LT_; // take a copy of LT_
const scalarSquareMatrix Avar = A_; // take a copy of A_
LT_ = scalarSquareMatrix(nActiveSpecies_+3, Zero);
A_ = scalarSquareMatrix(nActiveSpecies_+3, Zero);
LT_ = scalarSquareMatrix(nActive_+3, Zero);
A_ = scalarSquareMatrix(nActive_+3, Zero);
// write the initial active species
for (label i=0; i<initNActiveSpecies; i++)
@ -640,27 +637,27 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
{
for (label j=1; j>=0; j--)
{
LT_(i, nActiveSpecies_+j)=LTvar(i, initNActiveSpecies+j);
A_(i, nActiveSpecies_+j)=Avar(i, initNActiveSpecies+j);
LT_(nActiveSpecies_+j, i)=LTvar(initNActiveSpecies+j, i);
A_(nActiveSpecies_+j, i)=Avar(initNActiveSpecies+j, i);
LT_(i, nActive_+j)=LTvar(i, initNActiveSpecies+j);
A_(i, nActive_+j)=Avar(i, initNActiveSpecies+j);
LT_(nActive_+j, i)=LTvar(initNActiveSpecies+j, i);
A_(nActive_+j, i)=Avar(initNActiveSpecies+j, i);
}
}
// end with the diagonal elements for temperature and pressure
LT_(nActiveSpecies_, nActiveSpecies_)=
LT_(nActive_, nActive_)=
LTvar(initNActiveSpecies, initNActiveSpecies);
A_(nActiveSpecies_, nActiveSpecies_)=
A_(nActive_, nActive_)=
Avar(initNActiveSpecies, initNActiveSpecies);
LT_(nActiveSpecies_+1, nActiveSpecies_+1)=
LT_(nActive_+1, nActive_+1)=
LTvar(initNActiveSpecies+1, initNActiveSpecies+1);
A_(nActiveSpecies_+1, nActiveSpecies_+1)=
A_(nActive_+1, nActive_+1)=
Avar(initNActiveSpecies+1, initNActiveSpecies+1);
LT_(nActiveSpecies_+2, nActiveSpecies_+2)=
LT_(nActive_+2, nActive_+2)=
LTvar(initNActiveSpecies+2, initNActiveSpecies+2);
A_(nActiveSpecies_+2, nActiveSpecies_+2)=
A_(nActive_+2, nActive_+2)=
Avar(initNActiveSpecies+2, initNActiveSpecies+2);
for (label i=initNActiveSpecies; i<nActiveSpecies_;i++)
for (label i=initNActiveSpecies; i<nActive_;i++)
{
LT_(i, i)=
1.0
@ -671,7 +668,7 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
}
const label dim =
mechRedActive ? nActiveSpecies_ + 3 : completeSpaceSize();
table_.reduction() ? nActive_ + 3 : completeSpaceSize();
// beginning of grow algorithm
scalarField phiTilde(dim, 0);
@ -683,7 +680,7 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
for (label j=i; j<dim-3; j++)// LT is upper triangular
{
const label sj =
mechRedActive
table_.reduction()
? simplifiedToCompleteIndex_[j]
: j;
@ -747,19 +744,19 @@ Foam::label Foam::chemPointISAT<ThermoType>::simplifiedToCompleteIndex
const label i
)
{
if (i < nActiveSpecies_)
if (i < nActive_)
{
return simplifiedToCompleteIndex_[i];
}
else if (i == nActiveSpecies_)
else if (i == nActive_)
{
return completeSpaceSize_ - 3;
}
else if (i == nActiveSpecies_ + 1)
else if (i == nActive_ + 1)
{
return completeSpaceSize_ - 2;
}
else if (i == nActiveSpecies_ + 2)
else if (i == nActive_ + 2)
{
return completeSpaceSize_ - 1;
}

25
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -184,7 +184,7 @@ class chemPointISAT
static scalar tolerance_;
//- Number of active species stored in the chemPoint
label nActiveSpecies_;
label nActive_;
//- Vectors that store the index conversion between the simplified
// and the complete chemical mechanism
@ -256,8 +256,9 @@ public:
const scalarField& Rphi,
const scalarSquareMatrix& A,
const scalarField& scaleFactor,
const scalar& tolerance,
const label& completeSpaceSize,
const scalar tolerance,
const label completeSpaceSize,
const label nActive,
const dictionary& coeffsDict,
binaryNode<ThermoType>* node = nullptr
);
@ -275,13 +276,13 @@ public:
// Member Functions
//- Access to the chemistryModel
inline const chemistryModel<ThermoType>& chemistry()
//- Access to the ISAT table
inline const chemistryTabulationMethods::ISAT<ThermoType>& table() const
{
return table_.chemistry();
return table_;
}
inline label nGrowth()
inline label nGrowth() const
{
return nGrowth_;
}
@ -301,12 +302,12 @@ public:
return Rphi_;
}
inline const scalarField& scaleFactor()
inline const scalarField& scaleFactor() const
{
return scaleFactor_;
}
inline const scalar& tolerance()
inline const scalar& tolerance() const
{
return tolerance_;
}
@ -341,9 +342,9 @@ public:
return LT_;
}
inline label nActiveSpecies()
inline label nActive()
{
return nActiveSpecies_;
return nActive_;
}
inline List<label>& completeToSimplifiedIndex()

2
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H
View File

@ -56,7 +56,6 @@ class chemistryModel;
template<class ThermoType>
class chemistryTabulationMethod
{
public:
//- Runtime type information
@ -123,6 +122,7 @@ public:
(
const scalarField& phiQ,
const scalarField& RphiQ,
const label nActive,
const label li,
const scalar deltaT
) = 0;

1
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H
View File

@ -103,6 +103,7 @@ public:
(
const scalarField& phiq,
const scalarField& Rphiq,
const label nActive,
const label li,
const scalar deltaT
)