33 template<
class CompType,
class ThermoType>
38 template<
class CompType,
class ThermoType>
49 for (
label k=0; k<nCols-1; k++)
52 for (
label i=k; i<nCols; i++)
62 for (
label i=k; i<nCols; i++)
67 for (
label i=k; i<nCols; i++)
75 for (
label j=k+1; j<nCols; j++)
78 for (
label i=k; i<nCols; i++)
80 sum +=
R(i, k)*
R(i, j);
82 scalar tau = sum/c[
k];
83 for (
label i=k; i<nCols; i++)
85 R(i, j) -= tau*
R(i, k);
91 d[nCols-1] =
R(nCols-1, nCols-1);
94 for (
label i=0; i<nCols; i++)
97 for (
label j=0; j<i; j++)
105 template<
class CompType,
class ThermoType>
117 for (k=n-1; k>=0; k--)
130 for (
label i=k-1; i>=0; i--)
132 rotate(R, i, w[i],-w[i+1], n);
137 else if (
mag(w[i]) >
mag(w[i+1]))
139 w[i] =
mag(w[i])*
sqrt(1.0 +
sqr(w[i+1]/w[i]));
143 w[i] =
mag(w[i+1])*
sqrt(1.0 +
sqr(w[i]/w[i+1]));
149 R(0, i) += w[0]*v[i];
154 rotate(R, i,
R(i, i), -
R(i+1, i), n);
159 template<
class CompType,
class ThermoType>
169 scalar
c, fact,
s, w,
y;
173 s = (b >= 0 ? 1 : -1);
199 template<
class CompType,
class ThermoType>
207 const scalar& tolerance,
208 const label& completeSpaceSize,
213 chemistry_(chemistry),
217 scaleFactor_(scaleFactor),
219 completeSpaceSize_(completeSpaceSize),
221 nActiveSpecies_(chemistry.
mechRed()->NsSimp()),
222 simplifiedToCompleteIndex_(nActiveSpecies_),
223 timeTag_(chemistry_.timeSteps()),
224 lastTimeUsed_(chemistry_.timeSteps()),
230 completeToSimplifiedIndex_
232 completeSpaceSize - (2 + (variableTimeStep() == 1 ? 1 : 0))
235 tolerance_ = tolerance;
237 if (variableTimeStep())
239 nAdditionalEqns_ = 3;
240 iddeltaT_ = completeSpaceSize - 1;
241 scaleFactor_[iddeltaT_] *= phi_[iddeltaT_] / tolerance_;
245 nAdditionalEqns_ = 2;
246 iddeltaT_ = completeSpaceSize;
248 idT_ = completeSpaceSize - nAdditionalEqns_;
249 idp_ = completeSpaceSize - nAdditionalEqns_ + 1;
251 bool isMechRedActive = chemistry_.mechRed()->active();
254 for (
label i=0; i<completeSpaceSize-nAdditionalEqns_; i++)
256 completeToSimplifiedIndex_[i] =
259 for (
label i=0; i<nActiveSpecies_; i++)
261 simplifiedToCompleteIndex_[i] =
266 label reduOrCompDim = completeSpaceSize;
269 reduOrCompDim = nActiveSpecies_+nAdditionalEqns_;
279 for (
label i=0; i<reduOrCompDim; i++)
281 D[i] =
max(S[i], 0.5);
290 for (
label i=0; i<reduOrCompDim; i++)
292 for (
label j=0; j<reduOrCompDim; j++)
298 compi = simplifiedToCompleteIndex(i);
305 Atilde(i, j) /= (tolerance*scaleFactor[compi]);
314 qrDecompose(reduOrCompDim, LT_);
318 template<
class CompType,
class ThermoType>
345 if (variableTimeStep())
347 nAdditionalEqns_ = 3;
348 idT_ = completeSpaceSize() - 3;
349 idp_ = completeSpaceSize() - 2;
350 iddeltaT_ = completeSpaceSize() - 1;
354 nAdditionalEqns_ = 2;
355 idT_ = completeSpaceSize() - 2;
356 idp_ = completeSpaceSize() - 1;
357 iddeltaT_ = completeSpaceSize();
364 template<
class CompType,
class ThermoType>
368 bool isMechRedActive = chemistry_.mechRed()->active();
372 dim = nActiveSpecies_;
376 dim = completeSpaceSize() - nAdditionalEqns_;
382 for (
label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
392 ||(isMechRedActive && completeToSimplifiedIndex_[i] != -1)
395 label si=(isMechRedActive) ? completeToSimplifiedIndex_[i] : i;
397 for (
label j=si; j<dim; j++)
399 label sj=(isMechRedActive) ? simplifiedToCompleteIndex_[j] : j;
400 temp += LT_(si, j)*dphi[sj];
403 temp += LT_(si, dim)*dphi[idT_];
404 temp += LT_(si, dim+1)*dphi[idp_];
405 if (variableTimeStep())
407 temp += LT_(si, dim+2)*dphi[iddeltaT_];
412 temp = dphi[i]/(tolerance_*scaleFactor_[i]);
415 epsTemp +=
sqr(temp);
417 if (printProportion_)
424 if (variableTimeStep())
429 LT_(dim, dim)*dphi[idT_]
430 +LT_(dim, dim+1)*dphi[idp_]
431 +LT_(dim, dim+2)*dphi[iddeltaT_]
439 LT_(dim, dim)*dphi[idT_]
440 +LT_(dim, dim+1)*dphi[idp_]
445 if (variableTimeStep())
450 LT_(dim+1, dim+1)*dphi[idp_]
451 +LT_(dim+1, dim+2)*dphi[iddeltaT_]
456 epsTemp +=
sqr(LT_(dim+1, dim+1)*dphi[idp_]);
459 if (variableTimeStep())
461 epsTemp +=
sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
464 if (printProportion_)
468 LT_(dim, dim)*dphi[idT_]
469 + LT_(dim, dim+1)*dphi[idp_]
473 sqr(LT_(dim+1, dim+1)*dphi[idp_]);
475 if (variableTimeStep())
478 sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
483 if (
sqrt(epsTemp) > 1 + tolerance_)
485 if (printProportion_)
489 for (
label i=0; i<completeSpaceSize(); i++)
498 if (maxIndex >= completeSpaceSize() - nAdditionalEqns_)
500 if (maxIndex == idT_)
504 else if (maxIndex == idp_)
508 else if (maxIndex == iddeltaT_)
515 propName = chemistry_.Y()[maxIndex].name();
517 Info<<
"Direction maximum impact to error in ellipsoid: " 519 Info<<
"Proportion to the total error on the retrieve: " 520 << max/(epsTemp+SMALL) << endl;
531 template<
class CompType,
class ThermoType>
543 bool isMechRedActive = chemistry_.mechRed()->active();
545 label dim = completeSpaceSize()-2;
548 dim = nActiveSpecies_;
553 for (
label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
558 label si = completeToSimplifiedIndex_[i];
563 for (
label j=0; j<dim; j++)
565 label sj=simplifiedToCompleteIndex_[j];
566 dRl += Avar(si, j)*dphi[sj];
568 dRl += Avar(si, nActiveSpecies_)*dphi[idT_];
569 dRl += Avar(si, nActiveSpecies_+1)*dphi[idp_];
570 if (variableTimeStep())
572 dRl += Avar(si, nActiveSpecies_+2)*dphi[iddeltaT_];
582 for (
label j=0; j<completeSpaceSize(); j++)
584 dRl += Avar(i, j)*dphi[j];
587 eps2 +=
sqr((dR[i]-dRl)/scaleFactorV[i]);
591 if (eps2 > tolerance())
603 template<
class CompType,
class ThermoType>
607 label dim = completeSpaceSize();
608 label initNActiveSpecies(nActiveSpecies_);
609 bool isMechRedActive = chemistry_.mechRed()->active();
613 label activeAdded(0);
618 for (
label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
624 completeToSimplifiedIndex_[i] == -1
625 && chemistry_.completeToSimplifiedIndex()[i]!=-1
635 completeToSimplifiedIndex_[i] != -1
636 && chemistry_.completeToSimplifiedIndex()[i] == -1
647 completeToSimplifiedIndex_[i] == -1
648 && chemistry_.completeToSimplifiedIndex()[i] == -1
658 if (activeAdded > maxNumNewDim_)
664 nActiveSpecies_ += dimToAdd.
size();
665 simplifiedToCompleteIndex_.setSize(nActiveSpecies_);
668 label si = nActiveSpecies_ - dimToAdd.
size() + i;
670 simplifiedToCompleteIndex_[si] = dimToAdd[i];
671 completeToSimplifiedIndex_[dimToAdd[i]] = si;
680 if (nActiveSpecies_ > initNActiveSpecies)
688 for (
label i=0; i<initNActiveSpecies; i++)
690 for (
label j=0; j<initNActiveSpecies; j++)
692 LT_(i, j) = LTvar(i, j);
693 A_(i, j) = Avar(i, j);
698 for (
label i=0; i<initNActiveSpecies; i++)
700 for (
label j=1; j>=0; j--)
702 LT_(i, nActiveSpecies_+j)=LTvar(i, initNActiveSpecies+j);
703 A_(i, nActiveSpecies_+j)=Avar(i, initNActiveSpecies+j);
704 LT_(nActiveSpecies_+j, i)=LTvar(initNActiveSpecies+j, i);
705 A_(nActiveSpecies_+j, i)=Avar(initNActiveSpecies+j, i);
709 LT_(nActiveSpecies_, nActiveSpecies_)=
710 LTvar(initNActiveSpecies, initNActiveSpecies);
711 A_(nActiveSpecies_, nActiveSpecies_)=
712 Avar(initNActiveSpecies, initNActiveSpecies);
713 LT_(nActiveSpecies_+1, nActiveSpecies_+1)=
714 LTvar(initNActiveSpecies+1, initNActiveSpecies+1);
715 A_(nActiveSpecies_+1, nActiveSpecies_+1)=
716 Avar(initNActiveSpecies+1, initNActiveSpecies+1);
718 if (variableTimeStep())
720 LT_(nActiveSpecies_+2, nActiveSpecies_+2)=
721 LTvar(initNActiveSpecies+2, initNActiveSpecies+2);
722 A_(nActiveSpecies_+2, nActiveSpecies_+2)=
723 Avar(initNActiveSpecies+2, initNActiveSpecies+2);
726 for (
label i=initNActiveSpecies; i<nActiveSpecies_;i++)
730 /(tolerance_*scaleFactor_[simplifiedToCompleteIndex_[i]]);
735 dim = nActiveSpecies_ + nAdditionalEqns_;
740 scalar normPhiTilde = 0;
743 for (
label i=0; i<dim; i++)
745 for (
label j=i; j<dim-nAdditionalEqns_; j++)
750 sj = simplifiedToCompleteIndex_[j];
752 phiTilde[i] += LT_(i, j)*dphi[sj];
755 phiTilde[i] += LT_(i, dim-nAdditionalEqns_)*dphi[idT_];
756 phiTilde[i] += LT_(i, dim-nAdditionalEqns_+1)*dphi[idp_];
758 if (variableTimeStep())
760 phiTilde[i] += LT_(i, dim-nAdditionalEqns_ + 2)*dphi[iddeltaT_];
762 normPhiTilde +=
sqr(phiTilde[i]);
765 scalar invSqrNormPhiTilde = 1.0/normPhiTilde;
766 normPhiTilde =
sqrt(normPhiTilde);
769 scalar gamma = (1/normPhiTilde - 1)*invSqrNormPhiTilde;
773 for (
label i=0; i<dim; i++)
775 for (
label j=0; j<=i;j++)
777 v[i] += phiTilde[j]*LT_(j, i);
781 qrUpdate(LT_,dim, u, v);
788 template<
class CompType,
class ThermoType>
791 this->numRetrieve_++;
795 template<
class CompType,
class ThermoType>
798 this->numRetrieve_ = 0;
802 template<
class CompType,
class ThermoType>
809 template<
class CompType,
class ThermoType>
816 if (i < nActiveSpecies_)
818 return simplifiedToCompleteIndex_[i];
820 else if (i == nActiveSpecies_)
822 return completeSpaceSize_-nAdditionalEqns_;
824 else if (i == nActiveSpecies_ + 1)
826 return completeSpaceSize_-nAdditionalEqns_ + 1;
828 else if (variableTimeStep() && (i == nActiveSpecies_ + 2))
830 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)
Extends chemistryModel by adding the TDAC method.
#define forAll(list, i)
Loop across all elements in list.
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
const scalar & tolerance()
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const scalarSquareMatrix & A() const
const scalarSquareMatrix & LT() const
List< label > & simplifiedToCompleteIndex()
binaryNode< CompType, ThermoType > *& node()
A list of keyword definitions, which are a keyword followed by any number of values (e...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const scalarRectangularMatrix & U() const
Return U.
Form T() const
Return the transpose of the matrix.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
bool checkSolution(const scalarField &phiq, const scalarField &Rphiq)
If phiq is not in the EOA, then the mapping is computed.
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()
List< label > & completeToSimplifiedIndex()
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
label k
Boltzmann constant.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
DynamicList< label > & simplifiedToCompleteIndex()
const scalarDiagonalMatrix & S() const
Return the singular values.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
A class for handling words, derived from string.
const scalarField & scaleFactor()
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
bool grow(const scalarField &phiq)
More details about the minimum-volume ellipsoid covering an.
Singular value decomposition of a rectangular matrix.
const scalarField & phi() const
autoPtr< chemistryReductionMethod< CompType, ThermoType > > & mechRed()
volScalarField scalarField(fieldObject, mesh)
TDACChemistryModel< CompType, ThermoType > & chemistry()
Access to the TDACChemistryModel.
const scalarField & Rphi() const
const scalarRectangularMatrix & V() const
Return the square matrix V.
void resetNumRetrieve()
Resets the number of retrieves at each time step.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#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.
dimensioned< scalar > mag(const dimensioned< Type > &)
SquareMatrix< scalar > scalarSquareMatrix