chemPointISAT.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "chemPointISAT.H"
27 #include "SVD.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 // Defined as static to be able to dynamically change it during simulations
32 // (all chemPoints refer to the same object)
33 template<class CompType, class ThermoType>
35 
36 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
37 
38 template<class CompType, class ThermoType>
40 (
41  const label nCols,
43 )
44 {
45  scalarField c(nCols);
46  scalarField d(nCols);
47  scalar scale, sigma, sum;
48 
49  for (label k=0; k<nCols-1; k++)
50  {
51  scale = 0;
52  for (label i=k; i<nCols; i++)
53  {
54  scale=max(scale, mag(R(i, k)));
55  }
56  if (scale == 0)
57  {
58  c[k] = d[k] = 0;
59  }
60  else
61  {
62  for (label i=k; i<nCols; i++)
63  {
64  R(i, k) /= scale;
65  }
66  sum = 0;
67  for (label i=k; i<nCols; i++)
68  {
69  sum += sqr(R(i, k));
70  }
71  sigma = sign(R(k, k))*sqrt(sum);
72  R(k, k) += sigma;
73  c[k] = sigma*R(k, k);
74  d[k] = -scale*sigma;
75  for (label j=k+1; j<nCols; j++)
76  {
77  sum = 0;
78  for ( label i=k; i<nCols; i++)
79  {
80  sum += R(i, k)*R(i, j);
81  }
82  scalar tau = sum/c[k];
83  for ( label i=k; i<nCols; i++)
84  {
85  R(i, j) -= tau*R(i, k);
86  }
87  }
88  }
89  }
90 
91  d[nCols-1] = R(nCols-1, nCols-1);
92 
93  // form R
94  for (label i=0; i<nCols; i++)
95  {
96  R(i, i) = d[i];
97  for ( label j=0; j<i; j++)
98  {
99  R(i, j)=0;
100  }
101  }
102 }
103 
104 
105 template<class CompType, class ThermoType>
107 (
109  const label n,
110  const Foam::scalarField &u,
111  const Foam::scalarField &v
112 )
113 {
114  label k;
115 
116  scalarField w(u);
117  for (k=n-1; k>=0; k--)
118  {
119  if (w[k] != 0)
120  {
121  break;
122  }
123  }
124 
125  if (k < 0)
126  {
127  k = 0;
128  }
129 
130  for (label i=k-1; i>=0; i--)
131  {
132  rotate(R, i, w[i],-w[i+1], n);
133  if (w[i] == 0)
134  {
135  w[i] = mag(w[i+1]);
136  }
137  else if (mag(w[i]) > mag(w[i+1]))
138  {
139  w[i] = mag(w[i])*sqrt(1.0 + sqr(w[i+1]/w[i]));
140  }
141  else
142  {
143  w[i] = mag(w[i+1])*sqrt(1.0 + sqr(w[i]/w[i+1]));
144  }
145  }
146 
147  for (label i=0; i<n; i++)
148  {
149  R(0, i) += w[0]*v[i];
150  }
151 
152  for (label i=0; i<k; i++)
153  {
154  rotate(R, i, R(i, i), -R(i+1, i), n);
155  }
156 }
157 
158 
159 template<class CompType, class ThermoType>
161 (
163  const label i,
164  const scalar a,
165  const scalar b,
166  label n
167 )
168 {
169  scalar c, fact, s, w, y;
170  if (a == 0)
171  {
172  c = 0;
173  s = (b >= 0 ? 1 : -1);
174  }
175  else if (mag(a) > mag(b))
176  {
177  fact = b/a;
178  c = sign(a)/sqrt(1.0 + sqr(fact));
179  s = fact*c;
180  }
181  else
182  {
183  fact = a/b;
184  s = sign(b)/sqrt(1.0 + sqr(fact));
185  c = fact*s;
186  }
187  for (label j=i;j<n;j++)
188  {
189  y = R(i, j);
190  w = R(i+1, j);
191  R(i, j) = c*y-s*w;
192  R(i+1, j) = s*y+c*w;
193  }
194 }
195 
196 
197 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
198 
199 template<class CompType, class ThermoType>
201 (
203  const scalarField& phi,
204  const scalarField& Rphi,
205  const scalarSquareMatrix& A,
206  const scalarField& scaleFactor,
207  const scalar& tolerance,
208  const label& completeSpaceSize,
209  const dictionary& coeffsDict,
211 )
212 :
213  chemistry_(chemistry),
214  phi_(phi),
215  Rphi_(Rphi),
216  A_(A),
217  scaleFactor_(scaleFactor),
218  node_(node),
219  completeSpaceSize_(completeSpaceSize),
220  nGrowth_(0),
221  nActiveSpecies_(chemistry.mechRed()->NsSimp()),
222  simplifiedToCompleteIndex_(nActiveSpecies_),
223  timeTag_(chemistry_.timeSteps()),
224  lastTimeUsed_(chemistry_.timeSteps()),
225  toRemove_(false),
226  maxNumNewDim_(coeffsDict.lookupOrDefault("maxNumNewDim",0)),
227  printProportion_(coeffsDict.lookupOrDefault("printProportion",false)),
228  numRetrieve_(0),
229  nLifeTime_(0),
230  completeToSimplifiedIndex_
231  (
232  completeSpaceSize - (2 + (variableTimeStep() == 1 ? 1 : 0))
233  )
234 {
235  tolerance_ = tolerance;
236 
237  if (variableTimeStep())
238  {
239  nAdditionalEqns_ = 3;
240  iddeltaT_ = completeSpaceSize - 1;
241  scaleFactor_[iddeltaT_] *= phi_[iddeltaT_] / tolerance_;
242  }
243  else
244  {
245  nAdditionalEqns_ = 2;
246  iddeltaT_ = completeSpaceSize; // will not be used
247  }
248  idT_ = completeSpaceSize - nAdditionalEqns_;
249  idp_ = completeSpaceSize - nAdditionalEqns_ + 1;
250 
251  bool isMechRedActive = chemistry_.mechRed()->active();
252  if (isMechRedActive)
253  {
254  for (label i=0; i<completeSpaceSize-nAdditionalEqns_; i++)
255  {
256  completeToSimplifiedIndex_[i] =
257  chemistry.completeToSimplifiedIndex()[i];
258  }
259  for (label i=0; i<nActiveSpecies_; i++)
260  {
261  simplifiedToCompleteIndex_[i] =
262  chemistry.simplifiedToCompleteIndex()[i];
263  }
264  }
265 
266  label reduOrCompDim = completeSpaceSize;
267  if (isMechRedActive)
268  {
269  reduOrCompDim = nActiveSpecies_+nAdditionalEqns_;
270  }
271 
272  // SVD decomposition A = U*D*V^T
273  SVD svdA(A);
274 
275  scalarDiagonalMatrix D(reduOrCompDim);
276  const scalarDiagonalMatrix& S = svdA.S();
277 
278  // Replace the value of vector D by max(D, 1/2), first ISAT paper
279  for (label i=0; i<reduOrCompDim; i++)
280  {
281  D[i] = max(S[i], 0.5);
282  }
283 
284  // Rebuild A with max length, tol and scale factor before QR decomposition
285  scalarRectangularMatrix Atilde(reduOrCompDim);
286 
287  // Result stored in Atilde
288  multiply(Atilde, svdA.U(), D, svdA.V().T());
289 
290  for (label i=0; i<reduOrCompDim; i++)
291  {
292  for (label j=0; j<reduOrCompDim; j++)
293  {
294  label compi = i;
295 
296  if (isMechRedActive)
297  {
298  compi = simplifiedToCompleteIndex(i);
299  }
300 
301  // SF*A/tolerance
302  // (where SF is diagonal with inverse of scale factors)
303  // SF*A is the same as dividing each line by the scale factor
304  // corresponding to the species of this line
305  Atilde(i, j) /= (tolerance*scaleFactor[compi]);
306  }
307  }
308 
309  // The object LT_ (the transpose of the Q) describe the EOA, since we have
310  // A^T B^T B A that should be factorized into L Q^T Q L^T and is set in the
311  // qrDecompose function
312  LT_ = scalarSquareMatrix(Atilde);
313 
314  qrDecompose(reduOrCompDim, LT_);
315 }
316 
317 
318 template<class CompType, class ThermoType>
320 (
322 )
323 :
324  chemistry_(p.chemistry()),
325  phi_(p.phi()),
326  Rphi_(p.Rphi()),
327  LT_(p.LT()),
328  A_(p.A()),
329  scaleFactor_(p.scaleFactor()),
330  node_(p.node()),
331  completeSpaceSize_(p.completeSpaceSize()),
332  nGrowth_(p.nGrowth()),
333  nActiveSpecies_(p.nActiveSpecies()),
334  simplifiedToCompleteIndex_(p.simplifiedToCompleteIndex()),
335  timeTag_(p.timeTag()),
336  lastTimeUsed_(p.lastTimeUsed()),
337  toRemove_(p.toRemove()),
338  maxNumNewDim_(p.maxNumNewDim()),
339  numRetrieve_(0),
340  nLifeTime_(0),
341  completeToSimplifiedIndex_(p.completeToSimplifiedIndex())
342 {
343  tolerance_ = p.tolerance();
344 
345  if (variableTimeStep())
346  {
347  nAdditionalEqns_ = 3;
348  idT_ = completeSpaceSize() - 3;
349  idp_ = completeSpaceSize() - 2;
350  iddeltaT_ = completeSpaceSize() - 1;
351  }
352  else
353  {
354  nAdditionalEqns_ = 2;
355  idT_ = completeSpaceSize() - 2;
356  idp_ = completeSpaceSize() - 1;
357  iddeltaT_ = completeSpaceSize(); // will not be used
358  }
359 }
360 
361 
362 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
363 
364 template<class CompType, class ThermoType>
366 {
367  scalarField dphi(phiq-phi());
368  bool isMechRedActive = chemistry_.mechRed()->active();
369  label dim(0);
370  if (isMechRedActive)
371  {
372  dim = nActiveSpecies_;
373  }
374  else
375  {
376  dim = completeSpaceSize() - nAdditionalEqns_;
377  }
378 
379  scalar epsTemp = 0;
380  List<scalar> propEps(completeSpaceSize(), scalar(0));
381 
382  for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
383  {
384  scalar temp = 0;
385 
386  // When mechanism reduction is inactive OR on active species multiply L
387  // by dphi to get the distance in the active species direction else (for
388  // inactive species), just multiply the diagonal element and dphi
389  if
390  (
391  !(isMechRedActive)
392  ||(isMechRedActive && completeToSimplifiedIndex_[i] != -1)
393  )
394  {
395  label si=(isMechRedActive) ? completeToSimplifiedIndex_[i] : i;
396 
397  for (label j=si; j<dim; j++)// LT is upper triangular
398  {
399  label sj=(isMechRedActive) ? simplifiedToCompleteIndex_[j] : j;
400  temp += LT_(si, j)*dphi[sj];
401  }
402 
403  temp += LT_(si, dim)*dphi[idT_];
404  temp += LT_(si, dim+1)*dphi[idp_];
405  if (variableTimeStep())
406  {
407  temp += LT_(si, dim+2)*dphi[iddeltaT_];
408  }
409  }
410  else
411  {
412  temp = dphi[i]/(tolerance_*scaleFactor_[i]);
413  }
414 
415  epsTemp += sqr(temp);
416 
417  if (printProportion_)
418  {
419  propEps[i] = temp;
420  }
421  }
422 
423  // Temperature
424  if (variableTimeStep())
425  {
426  epsTemp +=
427  sqr
428  (
429  LT_(dim, dim)*dphi[idT_]
430  +LT_(dim, dim+1)*dphi[idp_]
431  +LT_(dim, dim+2)*dphi[iddeltaT_]
432  );
433  }
434  else
435  {
436  epsTemp +=
437  sqr
438  (
439  LT_(dim, dim)*dphi[idT_]
440  +LT_(dim, dim+1)*dphi[idp_]
441  );
442  }
443 
444  // Pressure
445  if (variableTimeStep())
446  {
447  epsTemp +=
448  sqr
449  (
450  LT_(dim+1, dim+1)*dphi[idp_]
451  +LT_(dim+1, dim+2)*dphi[iddeltaT_]
452  );
453  }
454  else
455  {
456  epsTemp += sqr(LT_(dim+1, dim+1)*dphi[idp_]);
457  }
458 
459  if (variableTimeStep())
460  {
461  epsTemp += sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
462  }
463 
464  if (printProportion_)
465  {
466  propEps[idT_] = sqr
467  (
468  LT_(dim, dim)*dphi[idT_]
469  + LT_(dim, dim+1)*dphi[idp_]
470  );
471 
472  propEps[idp_] =
473  sqr(LT_(dim+1, dim+1)*dphi[idp_]);
474 
475  if (variableTimeStep())
476  {
477  propEps[iddeltaT_] =
478  sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
479  }
480 
481  }
482 
483  if (sqrt(epsTemp) > 1 + tolerance_)
484  {
485  if (printProportion_)
486  {
487  scalar max = -1;
488  label maxIndex = -1;
489  for (label i=0; i<completeSpaceSize(); i++)
490  {
491  if(max < propEps[i])
492  {
493  max = propEps[i];
494  maxIndex = i;
495  }
496  }
497  word propName;
498  if (maxIndex >= completeSpaceSize() - nAdditionalEqns_)
499  {
500  if (maxIndex == idT_)
501  {
502  propName = "T";
503  }
504  else if (maxIndex == idp_)
505  {
506  propName = "p";
507  }
508  else if (maxIndex == iddeltaT_)
509  {
510  propName = "deltaT";
511  }
512  }
513  else
514  {
515  propName = chemistry_.Y()[maxIndex].member();
516  }
517  Info<< "Direction maximum impact to error in ellipsoid: "
518  << propName << endl;
519  Info<< "Proportion to the total error on the retrieve: "
520  << max/(epsTemp+small) << endl;
521  }
522  return false;
523  }
524  else
525  {
526  return true;
527  }
528 }
529 
530 
531 template<class CompType, class ThermoType>
533 (
534  const scalarField& phiq,
535  const scalarField& Rphiq
536 )
537 {
538  scalar eps2 = 0;
539  scalarField dR(Rphiq - Rphi());
540  scalarField dphi(phiq - phi());
541  const scalarField& scaleFactorV(scaleFactor());
542  const scalarSquareMatrix& Avar(A());
543  bool isMechRedActive = chemistry_.mechRed()->active();
544  scalar dRl = 0;
545  label dim = completeSpaceSize()-2;
546  if (isMechRedActive)
547  {
548  dim = nActiveSpecies_;
549  }
550 
551  // Since we build only the solution for the species, T and p are not
552  // included
553  for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
554  {
555  dRl = 0;
556  if (isMechRedActive)
557  {
558  label si = completeToSimplifiedIndex_[i];
559 
560  // If this species is active
561  if (si != -1)
562  {
563  for (label j=0; j<dim; j++)
564  {
565  label sj=simplifiedToCompleteIndex_[j];
566  dRl += Avar(si, j)*dphi[sj];
567  }
568  dRl += Avar(si, nActiveSpecies_)*dphi[idT_];
569  dRl += Avar(si, nActiveSpecies_+1)*dphi[idp_];
570  if (variableTimeStep())
571  {
572  dRl += Avar(si, nActiveSpecies_+2)*dphi[iddeltaT_];
573  }
574  }
575  else
576  {
577  dRl = dphi[i];
578  }
579  }
580  else
581  {
582  for (label j=0; j<completeSpaceSize(); j++)
583  {
584  dRl += Avar(i, j)*dphi[j];
585  }
586  }
587  eps2 += sqr((dR[i]-dRl)/scaleFactorV[i]);
588  }
589 
590  eps2 = sqrt(eps2);
591  if (eps2 > tolerance())
592  {
593  return false;
594  }
595  else
596  {
597  // if the solution is in the ellipsoid of accuracy
598  return true;
599  }
600 }
601 
602 
603 template<class CompType, class ThermoType>
605 {
606  scalarField dphi(phiq - phi());
607  label dim = completeSpaceSize();
608  label initNActiveSpecies(nActiveSpecies_);
609  bool isMechRedActive = chemistry_.mechRed()->active();
610 
611  if (isMechRedActive)
612  {
613  label activeAdded(0);
614  DynamicList<label> dimToAdd(0);
615 
616  // check if the difference of active species is lower than the maximum
617  // number of new dimensions allowed
618  for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
619  {
620  // first test if the current chemPoint has an inactive species
621  // corresponding to an active one in the query point
622  if
623  (
624  completeToSimplifiedIndex_[i] == -1
625  && chemistry_.completeToSimplifiedIndex()[i]!=-1
626  )
627  {
628  activeAdded++;
629  dimToAdd.append(i);
630  }
631  // then test if an active species in the current chemPoint
632  // corresponds to an inactive on the query side
633  if
634  (
635  completeToSimplifiedIndex_[i] != -1
636  && chemistry_.completeToSimplifiedIndex()[i] == -1
637  )
638  {
639  activeAdded++;
640  // we don't need to add a new dimension but we count it to have
641  // control on the difference through maxNumNewDim
642  }
643  // finally test if both points have inactive species but
644  // with a dphi!=0
645  if
646  (
647  completeToSimplifiedIndex_[i] == -1
648  && chemistry_.completeToSimplifiedIndex()[i] == -1
649  && dphi[i] != 0
650  )
651  {
652  activeAdded++;
653  dimToAdd.append(i);
654  }
655  }
656 
657  // if the number of added dimension is too large, growth fail
658  if (activeAdded > maxNumNewDim_)
659  {
660  return false;
661  }
662 
663  // the number of added dimension to the current chemPoint
664  nActiveSpecies_ += dimToAdd.size();
665  simplifiedToCompleteIndex_.setSize(nActiveSpecies_);
666  forAll(dimToAdd, i)
667  {
668  label si = nActiveSpecies_ - dimToAdd.size() + i;
669  // add the new active species
670  simplifiedToCompleteIndex_[si] = dimToAdd[i];
671  completeToSimplifiedIndex_[dimToAdd[i]] = si;
672  }
673 
674  // update LT and A :
675  //-add new column and line for the new active species
676  //-transfer last two lines of the previous matrix (p and T) to the end
677  // (change the diagonal position)
678  //-set all element of the new lines and columns to zero except diagonal
679  // (=1/(tolerance*scaleFactor))
680  if (nActiveSpecies_ > initNActiveSpecies)
681  {
682  scalarSquareMatrix LTvar = LT_; // take a copy of LT_
683  scalarSquareMatrix Avar = A_; // take a copy of A_
684  LT_ = scalarSquareMatrix(nActiveSpecies_+nAdditionalEqns_, Zero);
685  A_ = scalarSquareMatrix(nActiveSpecies_+nAdditionalEqns_, Zero);
686 
687  // write the initial active species
688  for (label i=0; i<initNActiveSpecies; i++)
689  {
690  for (label j=0; j<initNActiveSpecies; j++)
691  {
692  LT_(i, j) = LTvar(i, j);
693  A_(i, j) = Avar(i, j);
694  }
695  }
696 
697  // write the columns for temperature and pressure
698  for (label i=0; i<initNActiveSpecies; i++)
699  {
700  for (label j=1; j>=0; j--)
701  {
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);
706  }
707  }
708  // end with the diagonal elements for temperature and pressure
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);
717 
718  if (variableTimeStep())
719  {
720  LT_(nActiveSpecies_+2, nActiveSpecies_+2)=
721  LTvar(initNActiveSpecies+2, initNActiveSpecies+2);
722  A_(nActiveSpecies_+2, nActiveSpecies_+2)=
723  Avar(initNActiveSpecies+2, initNActiveSpecies+2);
724  }
725 
726  for (label i=initNActiveSpecies; i<nActiveSpecies_;i++)
727  {
728  LT_(i, i)=
729  1.0
730  /(tolerance_*scaleFactor_[simplifiedToCompleteIndex_[i]]);
731  A_(i, i) = 1;
732  }
733  }
734 
735  dim = nActiveSpecies_ + nAdditionalEqns_;
736  }
737 
738  // beginning of grow algorithm
739  scalarField phiTilde(dim, 0);
740  scalar normPhiTilde = 0;
741  // p' = L^T.(p-phi)
742 
743  for (label i=0; i<dim; i++)
744  {
745  for (label j=i; j<dim-nAdditionalEqns_; j++)// LT is upper triangular
746  {
747  label sj = j;
748  if (isMechRedActive)
749  {
750  sj = simplifiedToCompleteIndex_[j];
751  }
752  phiTilde[i] += LT_(i, j)*dphi[sj];
753  }
754 
755  phiTilde[i] += LT_(i, dim-nAdditionalEqns_)*dphi[idT_];
756  phiTilde[i] += LT_(i, dim-nAdditionalEqns_+1)*dphi[idp_];
757 
758  if (variableTimeStep())
759  {
760  phiTilde[i] += LT_(i, dim-nAdditionalEqns_ + 2)*dphi[iddeltaT_];
761  }
762  normPhiTilde += sqr(phiTilde[i]);
763  }
764 
765  scalar invSqrNormPhiTilde = 1.0/normPhiTilde;
766  normPhiTilde = sqrt(normPhiTilde);
767 
768  // gamma = (1/|p'| - 1)/|p'|^2
769  scalar gamma = (1/normPhiTilde - 1)*invSqrNormPhiTilde;
770  scalarField u(gamma*phiTilde);
771  scalarField v(dim, 0);
772 
773  for (label i=0; i<dim; i++)
774  {
775  for (label j=0; j<=i;j++)
776  {
777  v[i] += phiTilde[j]*LT_(j, i);
778  }
779  }
780 
781  qrUpdate(LT_,dim, u, v);
782  nGrowth_++;
783 
784  return true;
785 }
786 
787 
788 template<class CompType, class ThermoType>
790 {
791  this->numRetrieve_++;
792 }
793 
794 
795 template<class CompType, class ThermoType>
797 {
798  this->numRetrieve_ = 0;
799 }
800 
801 
802 template<class CompType, class ThermoType>
804 {
805  this->nLifeTime_++;
806 }
807 
808 
809 template<class CompType, class ThermoType>
812 (
813  const label i
814 )
815 {
816  if (i < nActiveSpecies_)
817  {
818  return simplifiedToCompleteIndex_[i];
819  }
820  else if (i == nActiveSpecies_)
821  {
822  return completeSpaceSize_-nAdditionalEqns_;
823  }
824  else if (i == nActiveSpecies_ + 1)
825  {
826  return completeSpaceSize_-nAdditionalEqns_ + 1;
827  }
828  else if (variableTimeStep() && (i == nActiveSpecies_ + 2))
829  {
830  return completeSpaceSize_-nAdditionalEqns_ + 2;
831  }
832  else
833  {
834  return -1;
835  }
836 }
837 
838 
839 // ************************************************************************* //
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)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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.
Definition: label.H:59
const scalarSquareMatrix & A() const
const scalarSquareMatrix & LT() const
surfaceScalarField & phi
List< label > & simplifiedToCompleteIndex()
binaryNode< CompType, ThermoType > *& node()
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const scalarRectangularMatrix & U() const
Return U.
Definition: SVDI.H:38
Form T() const
Return the transpose of the matrix.
Definition: Matrix.C:264
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
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.
Definition: Ostream.H:256
Leaf of the binary tree. The chemPoint stores the composition &#39;phi&#39;, the mapping of this composition ...
Field< label > & completeToSimplifiedIndex()
List< label > & completeToSimplifiedIndex()
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
label k
Boltzmann constant.
const label & timeTag()
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
DynamicList< label > & simplifiedToCompleteIndex()
scalar y
const scalarDiagonalMatrix & S() const
Return the singular values.
Definition: SVDI.H:50
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].
Definition: createFields.H:27
A class for handling words, derived from string.
Definition: word.H:59
const scalarField & scaleFactor()
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
bool grow(const scalarField &phiq)
More details about the minimum-volume ellipsoid covering an.
Singular value decomposition of a rectangular matrix.
Definition: SVD.H:51
static const zero Zero
Definition: zero.H:97
Node of the binary tree.
Definition: binaryNode.H:49
const scalarField & phi() const
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.
Definition: SVDI.H:44
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.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
SquareMatrix< scalar > scalarSquareMatrix