36template<
class CompType,
class ThermoType>
41template<
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)
109template<
class CompType,
class ThermoType>
112 scalarSquareMatrix&
R,
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);
163template<
class CompType,
class ThermoType>
166 scalarSquareMatrix&
R,
173 scalar
c, fact,
s, w,
y;
177 s = (
b >= 0 ? 1 : -1);
192 for (label j=i; j<
n ;++j)
204template<
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);
295 for (label i=0; i<reduOrCompDim; ++i)
297 for (label j=0; j<reduOrCompDim; ++j)
319 qrDecompose(reduOrCompDim, LT_);
323template<
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;
369template<
class CompType,
class ThermoType>
373 bool isMechRedActive = chemistry_.mechRed()->active();
377 dim = nActiveSpecies_;
381 dim = completeSpaceSize() - nAdditionalEqns_;
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;
533template<
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())
605template<
class CompType,
class ThermoType>
609 label dim = completeSpaceSize();
610 label initNActiveSpecies(nActiveSpecies_);
611 bool isMechRedActive = chemistry_.mechRed()->active();
615 label activeAdded(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);
791template<
class CompType,
class ThermoType>
794 this->numRetrieve_++;
798template<
class CompType,
class ThermoType>
801 this->numRetrieve_ = 0;
805template<
class CompType,
class ThermoType>
812template<
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;
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
#define R(A, B, C, D, E, F, K, M)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Form T() const
Return conjugate transpose of Matrix.
Singular value decomposition of a rectangular matrix.
const scalarRectangularMatrix & V() const
Return the square matrix V.
const scalarRectangularMatrix & U() const
Return U.
const scalarDiagonalMatrix & S() const
Return the singular values.
Extends StandardChemistryModel by adding the TDAC method.
autoPtr< chemistryReductionMethod< ReactionThermo, ThermoType > > & mechRed()
DynamicList< label > & simplifiedToCompleteIndex()
Field< label > & completeToSimplifiedIndex()
void size(const label n)
Older name for setAddressableSize.
Leaf of the binary tree. The chemPoint stores the composition 'phi', the mapping of this composition ...
void resetNumRetrieve()
Resets the number of retrieves at each time step.
bool inEOA(const scalarField &phiq)
To RETRIEVE the mapping from the stored chemPoint phi, the query.
void increaseNumRetrieve()
Increases the number of retrieves the chempoint has generated.
const scalarField & scaleFactor()
label & completeSpaceSize()
const scalar & tolerance()
bool grow(const scalarField &phiq)
More details about the minimum-volume ellipsoid covering an.
void increaseNLifeTime()
Increases the "counter" of the chP life.
bool checkSolution(const scalarField &phiq, const scalarField &Rphiq)
If phiq is not in the EOA, then the mapping is computed.
bool variableTimeStep() const
List< label > & simplifiedToCompleteIndex()
TDACChemistryModel< CompType, ThermoType > & chemistry()
Access to the TDACChemistryModel.
const scalarSquareMatrix & A() const
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
A class for handling words, derived from Foam::string.
BasicChemistryModel< psiReactionThermo > & chemistry
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))
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
const dimensionedScalar c
Speed of light in a vacuum.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
SquareMatrix< scalar > scalarSquareMatrix
static constexpr const zero Zero
Global zero (0)
constexpr char nl
The newline '\n' character (0x0a)
const dimensionedScalar & D
#define forAll(list, i)
Loop across all elements in list.