diff --git a/src/thermophysicalModels/chemistryModel/Make/options b/src/thermophysicalModels/chemistryModel/Make/options index de52d7e6ae..5997147f12 100644 --- a/src/thermophysicalModels/chemistryModel/Make/options +++ b/src/thermophysicalModels/chemistryModel/Make/options @@ -7,7 +7,8 @@ EXE_INC = \ -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ -I$(LIB_SRC)/turbulenceModels/compressible/lnInclude \ -I$(LIB_SRC)/ODE/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude LIB_LIBS = \ -lcompressibleTransportModels \ @@ -16,4 +17,5 @@ LIB_LIBS = \ -lspecie \ -lthermophysicalFunctions \ -lODE \ - -lfiniteVolume + -lfiniteVolume \ + -lmeshTools diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C index 21198d9d22..5fef714490 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C @@ -25,6 +25,7 @@ License #include "TDACChemistryModel.H" #include "UniformField.H" +#include "localEulerDdtScheme.H" #include "clockTime.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -37,6 +38,11 @@ Foam::TDACChemistryModel::TDACChemistryModel ) : chemistryModel(mesh, phaseName), + variableTimeStep_ + ( + mesh.time().controlDict().lookupOrDefault("adjustTimeStep", false) + || fv::localEulerDdt::enabled(mesh) + ), timeSteps_(0), NsDAC_(this->nSpecie_), completeC_(this->nSpecie_, 0), @@ -598,7 +604,7 @@ Foam::scalar Foam::TDACChemistryModel::solve const bool reduced = mechRed_->active(); - label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1:0); + label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0); basicMultiComponentMixture& composition = this->thermo().composition(); @@ -661,10 +667,10 @@ Foam::scalar Foam::TDACChemistryModel::solve phiq[i] = this->Y()[i][celli]; } phiq[this->nSpecie()]=Ti; - phiq[this->nSpecie()+1]=pi; + phiq[this->nSpecie() + 1]=pi; if (tabulation_->variableTimeStep()) { - phiq[this->nSpecie()+2] = deltaT[celli]; + phiq[this->nSpecie() + 2] = deltaT[celli]; } diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H index fd92c6a641..bf07574a3a 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H @@ -85,6 +85,8 @@ class TDACChemistryModel { // Private member data + bool variableTimeStep_; + label timeSteps_; // Mechanism reduction @@ -159,11 +161,11 @@ public: // Member Functions - inline label timeSteps() const - { - return timeSteps_; - } + //- Return true if the time-step is variable and/or non-uniform + inline bool variableTimeStep() const; + //- Return the number of chemistry evaluations (used by ISAT) + inline label timeSteps() const; //- Create and return a TDAC log file of the given name inline autoPtr logFile(const word& name) const; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H index 8aedea0571..362820308d 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,6 +25,22 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +inline bool +Foam::TDACChemistryModel::variableTimeStep() const +{ + return variableTimeStep_; +} + + +template +inline Foam::label +Foam::TDACChemistryModel::timeSteps() const +{ + return timeSteps_; +} + + template inline Foam::autoPtr Foam::TDACChemistryModel::logFile(const word& name) const diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C index 41de9be561..555365b7ba 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C @@ -98,30 +98,13 @@ Foam::chemistryTabulationMethods::ISAT::ISAT } } 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; @@ -433,7 +416,7 @@ void Foam::chemistryTabulationMethods::ISAT::computeA // For temperature and pressure, only unity on the diagonal 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; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.C index bfb092a737..5a4c4663b1 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.C @@ -28,19 +28,15 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -Foam::binaryNode::binaryNode -( -) +Foam::binaryNode::binaryNode() : leafLeft_(nullptr), leafRight_(nullptr), nodeLeft_(nullptr), nodeRight_(nullptr), parent_(nullptr), - variableTimeStep_(false), nAdditionalEqns_(0) -{ -} +{} template @@ -56,10 +52,9 @@ Foam::binaryNode::binaryNode nodeLeft_(nullptr), nodeRight_(nullptr), parent_(parent), - variableTimeStep_(elementLeft->variableTimeStep()), v_(elementLeft->completeSpaceSize(), 0) { - if (this->variableTimeStep_) + if (elementLeft->variableTimeStep()) { nAdditionalEqns_ = 3; } @@ -72,41 +67,11 @@ Foam::binaryNode::binaryNode a_ = calcA(elementLeft, elementRight); } -template -Foam::binaryNode::binaryNode -( - binaryNode *bn -) -: - leafLeft_(bn->leafLeft()), - leafRight_(bn->leafRight()), - nodeLeft_(bn->nodeLeft()), - nodeRight_(bn->nodeRight()), - parent_(bn->parent()), - variableTimeStep_ - ( - this->coeffsDict_.lookupOrDefault("variableTimeStep", false) - ), - v_(bn->v()), - a_(bn->a()) -{ - if (this->variableTimeStep_) - { - nAdditionalEqns_ = 3; - } - else - { - nAdditionalEqns_ = 2; - } -} - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - template -void -Foam::binaryNode::calcV +void Foam::binaryNode::calcV ( chemPointISAT*& elementLeft, chemPointISAT*& elementRight, @@ -116,7 +81,8 @@ Foam::binaryNode::calcV // LT is the transpose of the L matrix scalarSquareMatrix& LT = elementLeft->LT(); bool mechReductionActive = elementLeft->chemistry().mechRed()->active(); - // difference of composition in the full species domain + + // Difference of composition in the full species domain scalarField phiDif(elementRight->phi() - elementLeft->phi()); const scalarField& scaleFactor(elementLeft->scaleFactor()); scalar epsTol = elementLeft->tolerance(); @@ -131,9 +97,9 @@ Foam::binaryNode::calcV if (icompleteSpaceSize() - nAdditionalEqns_) { si = elementLeft->completeToSimplifiedIndex()[i]; - outOfIndexI = (si==-1); + outOfIndexI = (si == -1); } - else// temperature and pressure + else // temperature and pressure { outOfIndexI = false; const label dif = @@ -143,7 +109,7 @@ Foam::binaryNode::calcV } if (!mechReductionActive || (mechReductionActive && !(outOfIndexI))) { - v[i]=0.0; + v[i] = 0; for (label j=0; jcompleteSpaceSize(); j++) { label sj = j; @@ -173,7 +139,7 @@ Foam::binaryNode::calcV ||(mechReductionActive && !(outOfIndexJ)) ) { - // since L is a lower triangular matrix k=0->min(i, j) + // Since L is a lower triangular matrix k=0->min(i, j) for (label k=0; k<=min(si, sj); k++) { v[i] += LT(k, si)*LT(k, sj)*phiDif[j]; @@ -183,8 +149,8 @@ Foam::binaryNode::calcV } else { - // when it is an inactive species the diagonal element of LT is - // 1/(scaleFactor*epsTol) + // When it is an inactive species the diagonal element of LT is + // 1/(scaleFactor*epsTol) v[i] = phiDif[i]/sqr(scaleFactor[i]*epsTol); } } @@ -198,13 +164,13 @@ Foam::scalar Foam::binaryNode::calcA chemPointISAT* elementRight ) { - scalar a = 0.0; - scalarField phih((elementLeft->phi()+elementRight->phi())/2); - label completeSpaceSize = elementLeft->completeSpaceSize(); - for (label i=0; iphi() + elementRight->phi())/2); + scalar a = 0; + forAll(phih, i) { a += v_[i]*phih[i]; } + return a; } diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.H index 4b23e805d7..99952a51f6 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.H @@ -67,9 +67,6 @@ public: //- Parent node binaryNode* parent_; - //- Switch to allow variable time step (off by default) - bool variableTimeStep_; - //- Number of equations in addition to the species eqs. label nAdditionalEqns_; @@ -137,11 +134,6 @@ public: chemPointISAT* elementRight, binaryNode* parent ); - //- Construct from another binary node - binaryNode - ( - binaryNode *bn - ); // Member functions diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C index ecc212e1eb..33bf768813 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C @@ -227,18 +227,14 @@ Foam::chemPointISAT::chemPointISAT printProportion_(coeffsDict.lookupOrDefault("printProportion",false)), numRetrieve_(0), nLifeTime_(0), - variableTimeStep_ - ( - coeffsDict.lookupOrDefault("variableTimeStep", false) - ), completeToSimplifiedIndex_ ( - completeSpaceSize - (2 + (variableTimeStep_ == 1 ? 1 : 0)) + completeSpaceSize - (2 + (variableTimeStep() == 1 ? 1 : 0)) ) { - tolerance_=tolerance; + tolerance_ = tolerance; - if (this->variableTimeStep_) + if (variableTimeStep()) { nAdditionalEqns_ = 3; iddeltaT_ = completeSpaceSize - 1; @@ -342,12 +338,11 @@ Foam::chemPointISAT::chemPointISAT maxNumNewDim_(p.maxNumNewDim()), numRetrieve_(0), nLifeTime_(0), - variableTimeStep_(p.variableTimeStep()), completeToSimplifiedIndex_(p.completeToSimplifiedIndex()) { tolerance_ = p.tolerance(); - if (this->variableTimeStep_) + if (variableTimeStep()) { nAdditionalEqns_ = 3; idT_ = completeSpaceSize() - 3; @@ -407,7 +402,7 @@ bool Foam::chemPointISAT::inEOA(const scalarField& phiq) temp += LT_(si, dim)*dphi[idT_]; temp += LT_(si, dim+1)*dphi[idp_]; - if (variableTimeStep_) + if (variableTimeStep()) { temp += LT_(si, dim+2)*dphi[iddeltaT_]; } @@ -426,7 +421,7 @@ bool Foam::chemPointISAT::inEOA(const scalarField& phiq) } // Temperature - if (variableTimeStep_) + if (variableTimeStep()) { epsTemp += sqr @@ -447,7 +442,7 @@ bool Foam::chemPointISAT::inEOA(const scalarField& phiq) } // Pressure - if (variableTimeStep_) + if (variableTimeStep()) { epsTemp += sqr @@ -461,7 +456,7 @@ bool Foam::chemPointISAT::inEOA(const scalarField& phiq) epsTemp += sqr(LT_(dim+1, dim+1)*dphi[idp_]); } - if (variableTimeStep_) + if (variableTimeStep()) { epsTemp += sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]); } @@ -477,7 +472,7 @@ bool Foam::chemPointISAT::inEOA(const scalarField& phiq) propEps[idp_] = sqr(LT_(dim+1, dim+1)*dphi[idp_]); - if (variableTimeStep_) + if (variableTimeStep()) { propEps[iddeltaT_] = sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]); @@ -572,7 +567,7 @@ bool Foam::chemPointISAT::checkSolution } dRl += Avar(si, nActiveSpecies_)*dphi[idT_]; dRl += Avar(si, nActiveSpecies_+1)*dphi[idp_]; - if (variableTimeStep_) + if (variableTimeStep()) { dRl += Avar(si, nActiveSpecies_+2)*dphi[iddeltaT_]; } @@ -719,7 +714,8 @@ bool Foam::chemPointISAT::grow(const scalarField& phiq) LTvar(initNActiveSpecies+1, initNActiveSpecies+1); A_(nActiveSpecies_+1, nActiveSpecies_+1)= Avar(initNActiveSpecies+1, initNActiveSpecies+1); - if (variableTimeStep_) + + if (variableTimeStep()) { LT_(nActiveSpecies_+2, nActiveSpecies_+2)= LTvar(initNActiveSpecies+2, initNActiveSpecies+2); @@ -755,9 +751,11 @@ bool Foam::chemPointISAT::grow(const scalarField& phiq) } phiTilde[i] += LT_(i, j)*dphi[sj]; } + phiTilde[i] += LT_(i, dim-nAdditionalEqns_)*dphi[idT_]; phiTilde[i] += LT_(i, dim-nAdditionalEqns_+1)*dphi[idp_]; - if (variableTimeStep_) + + if (variableTimeStep()) { phiTilde[i] += LT_(i, dim-nAdditionalEqns_ + 2)*dphi[iddeltaT_]; } @@ -827,7 +825,7 @@ simplifiedToCompleteIndex { return completeSpaceSize_-nAdditionalEqns_ + 1; } - else if (variableTimeStep_ && (i == nActiveSpecies_ + 2)) + else if (variableTimeStep() && (i == nActiveSpecies_ + 2)) { return completeSpaceSize_-nAdditionalEqns_ + 2; } diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H index cae8ea213b..3b9ed78056 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H @@ -198,9 +198,6 @@ class chemPointISAT // to still live according to the maxChPLifeTime_ parameter label nLifeTime_; - //- Switch to allow variable time step (off by default) - Switch variableTimeStep_; - List