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))
240 tolerance_ = tolerance;
242 if (variableTimeStep())
244 nAdditionalEqns_ = 3;
245 iddeltaT_ = completeSpaceSize - 1;
246 scaleFactor_[iddeltaT_] *= phi_[iddeltaT_] / tolerance_;
250 nAdditionalEqns_ = 2;
251 iddeltaT_ = completeSpaceSize;
253 idT_ = completeSpaceSize - nAdditionalEqns_;
254 idp_ = completeSpaceSize - nAdditionalEqns_ + 1;
256 bool isMechRedActive = chemistry_.mechRed()->active();
259 for (label i=0; i<completeSpaceSize-nAdditionalEqns_; ++i)
261 completeToSimplifiedIndex_[i] =
262 chemistry.completeToSimplifiedIndex()[i];
264 for (label i=0; i<nActiveSpecies_; ++i)
266 simplifiedToCompleteIndex_[i] =
267 chemistry.simplifiedToCompleteIndex()[i];
271 label reduOrCompDim = completeSpaceSize;
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)
303 compi = simplifiedToCompleteIndex(i);
310 Atilde(i, j) /= (tolerance*scaleFactor[compi]);
319 qrDecompose(reduOrCompDim, LT_);
323 template<
class CompType,
class ThermoType>
329 chemistry_(
p.chemistry()),
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();
350 if (variableTimeStep())
352 nAdditionalEqns_ = 3;
353 idT_ = completeSpaceSize() - 3;
354 idp_ = completeSpaceSize() - 2;
355 iddeltaT_ = completeSpaceSize() - 1;
359 nAdditionalEqns_ = 2;
360 idT_ = completeSpaceSize() - 2;
361 idp_ = completeSpaceSize() - 1;
362 iddeltaT_ = completeSpaceSize();
369 template<
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;
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);
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;