36 template<
class CompType,
class ThermoType>
41 template<
class CompType,
class ThermoType>
52 for (label
k=0;
k<nCols-1; ++
k)
55 for (label i=
k; i<nCols; ++i)
66 for (label i=
k; i<nCols; ++i)
71 for (label i=
k; i<nCols; ++i)
79 for (label j=
k+1; j<nCols; ++j)
82 for ( label i=
k; i<nCols; ++i)
86 scalar tau =
sum/
c[
k];
87 for ( label i=
k; i<nCols; ++i)
89 R(i, j) -= tau*
R(i,
k);
95 d[nCols-1] =
R(nCols-1, nCols-1);
98 for (label i=0; i<nCols; ++i)
101 for ( label j=0; j<i; ++j)
109 template<
class CompType,
class ThermoType>
121 for (
k=
n-1;
k>=0; --
k)
134 for (label i=
k-1; i>=0; --i)
136 rotate(
R, i, w[i],-w[i+1],
n);
141 else if (
mag(w[i]) >
mag(w[i+1]))
143 w[i] =
mag(w[i])*
sqrt(1.0 +
sqr(w[i+1]/w[i]));
147 w[i] =
mag(w[i+1])*
sqrt(1.0 +
sqr(w[i]/w[i+1]));
151 for (label i=0; i<
n; ++i)
153 R(0, i) += w[0]*v[i];
156 for (label i=0; i<
k; ++i)
158 rotate(
R, i,
R(i, i), -
R(i+1, i),
n);
163 template<
class CompType,
class ThermoType>
173 scalar
c, fact,
s, w,
y;
177 s = (
b >= 0 ? 1 : -1);
192 for (label j=i; j<
n ;++j)
204 template<
class CompType,
class ThermoType>
212 const scalar& tolerance,
213 const label& completeSpaceSize,
222 scaleFactor_(scaleFactor),
224 completeSpaceSize_(completeSpaceSize),
226 nActiveSpecies_(
chemistry.mechRed()->NsSimp()),
227 simplifiedToCompleteIndex_(nActiveSpecies_),
228 timeTag_(chemistry_.timeSteps()),
229 lastTimeUsed_(chemistry_.timeSteps()),
231 maxNumNewDim_(coeffsDict.getOrDefault(
"maxNumNewDim", 0)),
232 printProportion_(coeffsDict.getOrDefault(
"printProportion", false)),
235 completeToSimplifiedIndex_
237 completeSpaceSize - (2 + (variableTimeStep() == 1 ? 1 : 0))
244 nAdditionalEqns_ = 3;
246 scaleFactor_[iddeltaT_] *= phi_[iddeltaT_] / tolerance_;
250 nAdditionalEqns_ = 2;
256 bool isMechRedActive = chemistry_.
mechRed()->active();
261 completeToSimplifiedIndex_[i] =
264 for (label i=0; i<nActiveSpecies_; ++i)
266 simplifiedToCompleteIndex_[i] =
274 reduOrCompDim = nActiveSpecies_+nAdditionalEqns_;
284 for (label i=0; i<reduOrCompDim; ++i)
286 D[i] =
max(S[i], 0.5);
293 multiply(Atilde, svdA.U(),
D, svdA.V().T());
295 for (label i=0; i<reduOrCompDim; ++i)
297 for (label j=0; j<reduOrCompDim; ++j)
319 qrDecompose(reduOrCompDim, LT_);
323 template<
class CompType,
class ThermoType>
334 scaleFactor_(
p.scaleFactor()),
336 completeSpaceSize_(
p.completeSpaceSize()),
337 nGrowth_(
p.nGrowth()),
338 nActiveSpecies_(
p.nActiveSpecies()),
339 simplifiedToCompleteIndex_(
p.simplifiedToCompleteIndex()),
340 timeTag_(
p.timeTag()),
341 lastTimeUsed_(
p.lastTimeUsed()),
342 toRemove_(
p.toRemove()),
343 maxNumNewDim_(
p.maxNumNewDim()),
346 completeToSimplifiedIndex_(
p.completeToSimplifiedIndex())
348 tolerance_ =
p.tolerance();
352 nAdditionalEqns_ = 3;
359 nAdditionalEqns_ = 2;
369 template<
class CompType,
class ThermoType>
373 bool isMechRedActive = chemistry_.mechRed()->active();
377 dim = nActiveSpecies_;
381 dim = completeSpaceSize() - nAdditionalEqns_;
385 List<scalar> propEps(completeSpaceSize(),
Zero);
387 for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; ++i)
397 ||(isMechRedActive && completeToSimplifiedIndex_[i] != -1)
400 label si = isMechRedActive ? completeToSimplifiedIndex_[i] : i;
402 for (label j=si; j<dim; ++j)
404 label sj = isMechRedActive ? simplifiedToCompleteIndex_[j] : j;
405 temp += LT_(si, j)*dphi[sj];
408 temp += LT_(si, dim)*dphi[idT_];
409 temp += LT_(si, dim+1)*dphi[idp_];
410 if (variableTimeStep())
412 temp += LT_(si, dim+2)*dphi[iddeltaT_];
417 temp = dphi[i]/(tolerance_*scaleFactor_[i]);
420 epsTemp +=
sqr(temp);
422 if (printProportion_)
429 if (variableTimeStep())
434 LT_(dim, dim)*dphi[idT_]
435 +LT_(dim, dim+1)*dphi[idp_]
436 +LT_(dim, dim+2)*dphi[iddeltaT_]
444 LT_(dim, dim)*dphi[idT_]
445 +LT_(dim, dim+1)*dphi[idp_]
450 if (variableTimeStep())
455 LT_(dim+1, dim+1)*dphi[idp_]
456 +LT_(dim+1, dim+2)*dphi[iddeltaT_]
461 epsTemp +=
sqr(LT_(dim+1, dim+1)*dphi[idp_]);
464 if (variableTimeStep())
466 epsTemp +=
sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
469 if (printProportion_)
473 LT_(dim, dim)*dphi[idT_]
474 + LT_(dim, dim+1)*dphi[idp_]
477 propEps[idp_] =
sqr(LT_(dim+1, dim+1)*dphi[idp_]);
479 if (variableTimeStep())
481 propEps[iddeltaT_] =
sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
486 if (
sqrt(epsTemp) > 1 + tolerance_)
488 if (printProportion_)
492 for (label i=0; i<completeSpaceSize(); ++i)
494 if (
max < propEps[i])
501 if (maxIndex >= completeSpaceSize() - nAdditionalEqns_)
503 if (maxIndex == idT_)
507 else if (maxIndex == idp_)
511 else if (maxIndex == iddeltaT_)
518 propName = chemistry_.Y()[maxIndex].member();
521 Info<<
"Direction maximum impact to error in ellipsoid: " 523 <<
"Proportion to the total error on the retrieve: " 524 <<
max/(epsTemp+SMALL) <<
endl;
533 template<
class CompType,
class ThermoType>
545 bool isMechRedActive = chemistry_.mechRed()->active();
547 label dim = completeSpaceSize() - 2;
550 dim = nActiveSpecies_;
555 for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; ++i)
560 label si = completeToSimplifiedIndex_[i];
565 for (label j=0; j<dim; ++j)
567 label sj=simplifiedToCompleteIndex_[j];
568 dRl += Avar(si, j)*dphi[sj];
570 dRl += Avar(si, nActiveSpecies_)*dphi[idT_];
571 dRl += Avar(si, nActiveSpecies_+1)*dphi[idp_];
572 if (variableTimeStep())
574 dRl += Avar(si, nActiveSpecies_+2)*dphi[iddeltaT_];
584 for (label j=0; j<completeSpaceSize(); ++j)
586 dRl += Avar(i, j)*dphi[j];
589 eps2 +=
sqr((dR[i]-dRl)/scaleFactorV[i]);
593 if (eps2 > tolerance())
605 template<
class CompType,
class ThermoType>
609 label dim = completeSpaceSize();
610 label initNActiveSpecies(nActiveSpecies_);
611 bool isMechRedActive = chemistry_.mechRed()->active();
615 label activeAdded(0);
616 DynamicList<label> dimToAdd(0);
620 for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; ++i)
626 completeToSimplifiedIndex_[i] == -1
627 && chemistry_.completeToSimplifiedIndex()[i]!=-1
637 completeToSimplifiedIndex_[i] != -1
638 && chemistry_.completeToSimplifiedIndex()[i] == -1
649 completeToSimplifiedIndex_[i] == -1
650 && chemistry_.completeToSimplifiedIndex()[i] == -1
660 if (activeAdded > maxNumNewDim_)
666 const label nDimAdd = dimToAdd.size();
667 nActiveSpecies_ += nDimAdd;
668 simplifiedToCompleteIndex_.setSize(nActiveSpecies_);
671 label si = nActiveSpecies_ - nDimAdd + i;
673 simplifiedToCompleteIndex_[si] = dimToAdd[i];
674 completeToSimplifiedIndex_[dimToAdd[i]] = si;
683 if (nActiveSpecies_ > initNActiveSpecies)
691 for (label i=0; i<initNActiveSpecies; ++i)
693 for (label j=0; j<initNActiveSpecies; ++j)
695 LT_(i, j) = LTvar(i, j);
696 A_(i, j) = Avar(i, j);
701 for (label i=0; i<initNActiveSpecies; ++i)
703 for (label j=1; j>=0; --j)
705 LT_(i, nActiveSpecies_+j) = LTvar(i, initNActiveSpecies+j);
706 A_(i, nActiveSpecies_+j) = Avar(i, initNActiveSpecies+j);
707 LT_(nActiveSpecies_+j, i) = LTvar(initNActiveSpecies+j, i);
708 A_(nActiveSpecies_+j, i) = Avar(initNActiveSpecies+j, i);
712 LT_(nActiveSpecies_, nActiveSpecies_) =
713 LTvar(initNActiveSpecies, initNActiveSpecies);
714 A_(nActiveSpecies_, nActiveSpecies_) =
715 Avar(initNActiveSpecies, initNActiveSpecies);
716 LT_(nActiveSpecies_+1, nActiveSpecies_+1) =
717 LTvar(initNActiveSpecies+1, initNActiveSpecies+1);
718 A_(nActiveSpecies_+1, nActiveSpecies_+1) =
719 Avar(initNActiveSpecies+1, initNActiveSpecies+1);
721 if (variableTimeStep())
723 LT_(nActiveSpecies_+2, nActiveSpecies_+2) =
724 LTvar(initNActiveSpecies+2, initNActiveSpecies+2);
725 A_(nActiveSpecies_+2, nActiveSpecies_+2) =
726 Avar(initNActiveSpecies+2, initNActiveSpecies+2);
729 for (label i=initNActiveSpecies; i<nActiveSpecies_; ++i)
733 /(tolerance_*scaleFactor_[simplifiedToCompleteIndex_[i]]);
738 dim = nActiveSpecies_ + nAdditionalEqns_;
743 scalar normPhiTilde = 0;
746 for (label i=0; i<dim; ++i)
748 for (label j=i; j<dim-nAdditionalEqns_; ++j)
753 sj = simplifiedToCompleteIndex_[j];
755 phiTilde[i] += LT_(i, j)*dphi[sj];
758 phiTilde[i] += LT_(i, dim-nAdditionalEqns_)*dphi[idT_];
759 phiTilde[i] += LT_(i, dim-nAdditionalEqns_+1)*dphi[idp_];
761 if (variableTimeStep())
763 phiTilde[i] += LT_(i, dim-nAdditionalEqns_ + 2)*dphi[iddeltaT_];
765 normPhiTilde +=
sqr(phiTilde[i]);
768 scalar invSqrNormPhiTilde = 1.0/normPhiTilde;
769 normPhiTilde =
sqrt(normPhiTilde);
772 scalar
gamma = (1/normPhiTilde - 1)*invSqrNormPhiTilde;
776 for (label i=0; i<dim; ++i)
778 for (label j=0; j<=i; ++j)
780 v[i] += phiTilde[j]*LT_(j, i);
784 qrUpdate(LT_, dim, u, v);
791 template<
class CompType,
class ThermoType>
794 this->numRetrieve_++;
798 template<
class CompType,
class ThermoType>
801 this->numRetrieve_ = 0;
805 template<
class CompType,
class ThermoType>
812 template<
class CompType,
class ThermoType>
819 if (i < nActiveSpecies_)
821 return simplifiedToCompleteIndex_[i];
823 else if (i == nActiveSpecies_)
825 return completeSpaceSize_-nAdditionalEqns_;
827 else if (i == nActiveSpecies_ + 1)
829 return completeSpaceSize_-nAdditionalEqns_ + 1;
831 else if (variableTimeStep() && (i == nActiveSpecies_ + 2))
833 return completeSpaceSize_-nAdditionalEqns_ + 2;
chemPointISAT(TDACChemistryModel< CompType, ThermoType > &chemistry, const scalarField &phi, const scalarField &Rphi, const scalarSquareMatrix &A, const scalarField &scaleFactor, const scalar &tolerance, const label &completeSpaceSize, const dictionary &coeffsDict, binaryNode< CompType, ThermoType > *node=nullptr)
Construct from components.
dimensionedScalar sign(const dimensionedScalar &ds)
DiagonalMatrix< scalar > scalarDiagonalMatrix
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
const scalar & tolerance()
const scalarSquareMatrix & A() const
List< label > & simplifiedToCompleteIndex()
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A list of keyword definitions, which are a keyword followed by a number of values (eg...
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
bool checkSolution(const scalarField &phiq, const scalarField &Rphiq)
If phiq is not in the EOA, then the mapping is computed.
autoPtr< chemistryReductionMethod< ReactionThermo, ThermoType > > & mechRed()
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Leaf of the binary tree. The chemPoint stores the composition 'phi', the mapping of this composition ...
Field< label > & completeToSimplifiedIndex()
label k
Boltzmann constant.
BasicChemistryModel< psiReactionThermo > & chemistry
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
#define forAll(list, i)
Loop across all elements in list.
DynamicList< label > & simplifiedToCompleteIndex()
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
const scalarField & scaleFactor()
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
bool grow(const scalarField &phiq)
More details about the minimum-volume ellipsoid covering an.
bool variableTimeStep() const
TDACChemistryModel< CompType, ThermoType > & chemistry()
Access to the TDACChemistryModel.
void resetNumRetrieve()
Resets the number of retrieves at each time step.
#define R(A, B, C, D, E, F, K, M)
void increaseNLifeTime()
Increases the "counter" of the chP life.
label & completeSpaceSize()
const dimensionedScalar c
Speed of light in a vacuum.
void increaseNumRetrieve()
Increases the number of retrieves the chempoint has generated.
bool inEOA(const scalarField &phiq)
To RETRIEVE the mapping from the stored chemPoint phi, the query.
messageStream Info
Information stream (stdout output on master, null elsewhere)
const dimensionedScalar & D
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
SquareMatrix< scalar > scalarSquareMatrix
RectangularMatrix< scalar > scalarRectangularMatrix
static constexpr const zero Zero
Global zero (0)