All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ISAT.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-2020 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 "ISAT.H"
27 #include "LUscalarMatrix.H"
28 
29 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
30 
31 template<class CompType, class ThermoType>
33 (
34  const dictionary& chemistryProperties,
36 )
37 :
39  (
40  chemistryProperties,
41  chemistry
42  ),
43  chemisTree_(chemistry, this->coeffsDict_),
44  scaleFactor_(chemistry.nEqns() + ((this->variableTimeStep()) ? 1 : 0), 1),
45  runTime_(chemistry.time()),
46  chPMaxLifeTime_
47  (
48  this->coeffsDict_.lookupOrDefault("chPMaxLifeTime", INT_MAX)
49  ),
50  maxGrowth_(this->coeffsDict_.lookupOrDefault("maxGrowth", INT_MAX)),
51  checkEntireTreeInterval_
52  (
53  this->coeffsDict_.lookupOrDefault("checkEntireTreeInterval", INT_MAX)
54  ),
55  maxDepthFactor_
56  (
57  this->coeffsDict_.lookupOrDefault
58  (
59  "maxDepthFactor",
60  (chemisTree_.maxNLeafs() - 1)
61  /(log(scalar(chemisTree_.maxNLeafs()))/log(2.0))
62  )
63  ),
64  minBalanceThreshold_
65  (
66  this->coeffsDict_.lookupOrDefault
67  (
68  "minBalanceThreshold", 0.1*chemisTree_.maxNLeafs()
69  )
70  ),
71  MRURetrieve_(this->coeffsDict_.lookupOrDefault("MRURetrieve", false)),
72  maxMRUSize_(this->coeffsDict_.lookupOrDefault("maxMRUSize", 0)),
73  lastSearch_(nullptr),
74  growPoints_(this->coeffsDict_.lookupOrDefault("growPoints", true)),
75  nRetrieved_(0),
76  nGrowth_(0),
77  nAdd_(0),
78  cleaningRequired_(false)
79 {
80  if (this->active_)
81  {
82  dictionary scaleDict(this->coeffsDict_.subDict("scaleFactor"));
83  label Ysize = this->chemistry_.Y().size();
84  scalar otherScaleFactor = scaleDict.lookup<scalar>("otherSpecies");
85  for (label i=0; i<Ysize; i++)
86  {
87  if (!scaleDict.found(this->chemistry_.Y()[i].member()))
88  {
89  scaleFactor_[i] = otherScaleFactor;
90  }
91  else
92  {
93  scaleFactor_[i] =
94  scaleDict.lookup<scalar>(this->chemistry_.Y()[i].member());
95  }
96  }
97  scaleFactor_[Ysize] = scaleDict.lookup<scalar>("Temperature");
98  scaleFactor_[Ysize + 1] = scaleDict.lookup<scalar>("Pressure");
99  if (this->variableTimeStep())
100  {
101  scaleFactor_[Ysize + 2] = scaleDict.lookup<scalar>("deltaT");
102  }
103  }
104 
105  if (this->variableTimeStep())
106  {
107  nAdditionalEqns_ = 3;
108  }
109  else
110  {
111  nAdditionalEqns_ = 2;
112  }
113 
114  if (this->log())
115  {
116  nRetrievedFile_ = chemistry.logFile("found_isat.out");
117  nGrowthFile_ = chemistry.logFile("growth_isat.out");
118  nAddFile_ = chemistry.logFile("add_isat.out");
119  sizeFile_ = chemistry.logFile("size_isat.out");
120  }
121 }
122 
123 
124 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
125 
126 template<class CompType, class ThermoType>
128 {}
129 
130 
131 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
132 
133 template<class CompType, class ThermoType>
135 (
137 )
138 {
139  if (maxMRUSize_ > 0 && MRURetrieve_)
140  {
141  // First search if the chemPoint is already in the list
142  bool isInList = false;
143  typename SLList <chemPointISAT<CompType, ThermoType>*>::iterator iter =
144  MRUList_.begin();
145  for ( ; iter != MRUList_.end(); ++iter)
146  {
147  if (iter() == phi0)
148  {
149  isInList = true;
150  break;
151  }
152  }
153  // If it is in the list, then move it to front
154  if (isInList)
155  {
156  if (iter() != MRUList_.first())
157  {
158  // iter hold the position of the element to move
159  MRUList_.remove(iter);
160 
161  // Insert the element in front of the list
162  MRUList_.insert(phi0);
163  }
164  }
165  else // chemPoint not yet in the list, iter is last
166  {
167  if (MRUList_.size() == maxMRUSize_)
168  {
169  MRUList_.remove(iter);
170  MRUList_.insert(phi0);
171  }
172  else
173  {
174  MRUList_.insert(phi0);
175  }
176  }
177  }
178 }
179 
180 
181 template<class CompType, class ThermoType>
183 (
185  const scalarField& phiq,
186  scalarField& Rphiq
187 )
188 {
189  label nEqns = this->chemistry_.nEqns(); // Species, T, p
190  bool mechRedActive = this->chemistry_.mechRed()->active();
191  Rphiq = phi0->Rphi();
192  scalarField dphi(phiq-phi0->phi());
193  const scalarSquareMatrix& gradientsMatrix = phi0->A();
194  List<label>& completeToSimplified(phi0->completeToSimplifiedIndex());
195 
196  // Rphiq[i]=Rphi0[i]+A(i, j)dphi[j]
197  // where Aij is dRi/dphi_j
198  for (label i=0; i<nEqns-nAdditionalEqns_; i++)
199  {
200  if (mechRedActive)
201  {
202  label si = completeToSimplified[i];
203  // The species is active
204  if (si != -1)
205  {
206  for (label j=0; j<nEqns-2; j++)
207  {
208  label sj = completeToSimplified[j];
209  if (sj != -1)
210  {
211  Rphiq[i] += gradientsMatrix(si, sj)*dphi[j];
212  }
213  }
214  Rphiq[i] +=
215  gradientsMatrix(si, phi0->nActiveSpecies())*dphi[nEqns - 2];
216  Rphiq[i] +=
217  gradientsMatrix(si, phi0->nActiveSpecies() + 1)
218  *dphi[nEqns - 1];
219 
220  if (this->variableTimeStep())
221  {
222  Rphiq[i] +=
223  gradientsMatrix(si, phi0->nActiveSpecies() + 2)
224  *dphi[nEqns];
225  }
226 
227  // As we use an approximation of A, Rphiq should be checked for
228  // negative values
229  Rphiq[i] = max(0, Rphiq[i]);
230  }
231  // The species is not active A(i, j) = I(i, j)
232  else
233  {
234  Rphiq[i] += dphi[i];
235  Rphiq[i] = max(0, Rphiq[i]);
236  }
237  }
238  else // Mechanism reduction is not active
239  {
240  for (label j=0; j<nEqns; j++)
241  {
242  Rphiq[i] += gradientsMatrix(i, j)*dphi[j];
243  }
244  // As we use a first order gradient matrix, Rphiq should be checked
245  // for negative values
246  Rphiq[i] = max(0, Rphiq[i]);
247  }
248  }
249 }
250 
251 
252 template<class CompType, class ThermoType>
254 (
256  const scalarField& phiq,
257  const scalarField& Rphiq
258 )
259 {
260  // If the pointer to the chemPoint is nullptr, the function stops
261  if (!phi0)
262  {
263  return false;
264  }
265 
266  // Raise a flag when the chemPoint used has been grown more than the
267  // allowed number of time
268  if (phi0->nGrowth() > maxGrowth_)
269  {
270  cleaningRequired_ = true;
271  phi0->toRemove() = true;
272  return false;
273  }
274 
275  // If the solution RphiQ is still within the tolerance we try to grow it
276  // in some cases this might result in a failure and the grow function of
277  // the chemPoint returns false
278  if (phi0->checkSolution(phiq, Rphiq))
279  {
280  return phi0->grow(phiq);
281  }
282  // The actual solution and the approximation given by ISAT are too different
283  else
284  {
285  return false;
286  }
287 }
288 
289 
290 template<class CompType, class ThermoType>
291 bool
293 {
294  bool treeModified(false);
295 
296  // Check all chemPoints to see if we need to delete some of the chemPoints
297  // according to the elapsed time and number of growths
298  chemPointISAT<CompType, ThermoType>* x = chemisTree_.treeMin();
299  while(x != nullptr)
300  {
302  chemisTree_.treeSuccessor(x);
303 
304  scalar elapsedTimeSteps = this->chemistry_.timeSteps() - x->timeTag();
305 
306  if ((elapsedTimeSteps > chPMaxLifeTime_) || (x->nGrowth() > maxGrowth_))
307  {
308  chemisTree_.deleteLeaf(x);
309  treeModified = true;
310  }
311  x = xtmp;
312  }
313 
314  MRUList_.clear();
315 
316  // Check if the tree should be balanced according to criterion:
317  // -the depth of the tree bigger than a*log2(size), log2(size) being the
318  // ideal depth (e.g. 4 leafs can be stored in a tree of depth 2)
319  if
320  (
321  chemisTree_.size() > minBalanceThreshold_
322  && chemisTree_.depth() >
323  maxDepthFactor_*log(scalar(chemisTree_.size()))/log(2.0)
324  )
325  {
326  chemisTree_.balance();
327  treeModified = true;
328  }
329 
330  // Return a bool to specify if the tree structure has been modified and is
331  // now below the user specified limit (true if not full)
332  return (treeModified && !chemisTree_.isFull());
333 }
334 
335 
336 template<class CompType, class ThermoType>
338 (
340  const scalarField& Rphiq,
341  const label li,
342  const scalar rhoi,
343  const scalar dt
344 )
345 {
346  bool mechRedActive = this->chemistry_.mechRed()->active();
347  label speciesNumber = this->chemistry_.nSpecie();
348  scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
349  for (label i=0; i<speciesNumber; i++)
350  {
351  label s2c = i;
352  if (mechRedActive)
353  {
354  s2c = this->chemistry_.simplifiedToCompleteIndex()[i];
355  }
356  Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermos()[s2c].W();
357  }
358  Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
359  Rcq[speciesNumber + 1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
360  if (this->variableTimeStep())
361  {
362  Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
363  }
364 
365  // Aaa is computed implicitly,
366  // A is given by A = C(psi0, t0+dt), where C is obtained through solving
367  // d/dt C(psi0, t) = J(psi(t))C(psi0, t)
368  // If we solve it implicitly:
369  // (C(psi0, t0+dt) - C(psi0, t0))/dt = J(psi(t0+dt))C(psi0, t0+dt)
370  // The Jacobian is thus computed according to the mapping
371  // C(psi0,t0+dt)*(I-dt*J(psi(t0+dt))) = C(psi0, t0)
372  // A = C(psi0,t0)/(I-dt*J(psi(t0+dt)))
373  // where C(psi0,t0) = I
374  scalarField dcdt(speciesNumber + 2, Zero);
375  this->chemistry_.jacobian(runTime_.value(), Rcq, li, dcdt, A);
376 
377  // The jacobian is computed according to the molar concentration
378  // the following conversion allows the code to use A with mass fraction
379  for (label i=0; i<speciesNumber; i++)
380  {
381  label si = i;
382 
383  if (mechRedActive)
384  {
385  si = this->chemistry_.simplifiedToCompleteIndex()[i];
386  }
387 
388  for (label j=0; j<speciesNumber; j++)
389  {
390  label sj = j;
391  if (mechRedActive)
392  {
393  sj = this->chemistry_.simplifiedToCompleteIndex()[j];
394  }
395  A(i, j) *=
396  -dt*this->chemistry_.specieThermos()[si].W()
397  /this->chemistry_.specieThermos()[sj].W();
398  }
399 
400  A(i, i) += 1;
401  // Columns for pressure and temperature
402  A(i, speciesNumber) *=
403  -dt*this->chemistry_.specieThermos()[si].W()/rhoi;
404  A(i, speciesNumber + 1) *=
405  -dt*this->chemistry_.specieThermos()[si].W()/rhoi;
406  }
407 
408  // For the temperature and pressure lines, ddc(dTdt)
409  // should be converted in ddY(dTdt)
410  for (label i=0; i<speciesNumber; i++)
411  {
412  label si = i;
413  if (mechRedActive)
414  {
415  si = this->chemistry_.simplifiedToCompleteIndex()[i];
416  }
417 
418  A(speciesNumber, i) *=
419  -dt*rhoi/this->chemistry_.specieThermos()[si].W();
420  A(speciesNumber + 1, i) *=
421  -dt*rhoi/this->chemistry_.specieThermos()[si].W();
422  }
423 
424  A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
425 
426  A(speciesNumber + 1, speciesNumber + 1) =
427  -dt*A(speciesNumber + 1, speciesNumber + 1) + 1;
428 
429  if (this->variableTimeStep())
430  {
431  A(speciesNumber + 2, speciesNumber + 2) = 1;
432  }
433 
434  // Inverse of (I-dt*J(psi(t0+dt)))
435  LUscalarMatrix LUA(A);
436  LUA.inv(A);
437 
438  // After inversion, lines of p and T are set to 0 except diagonal. This
439  // avoid skewness of the ellipsoid of accuracy and potential issues in the
440  // binary tree.
441  for (label i=0; i<speciesNumber; i++)
442  {
443  A(speciesNumber, i) = 0;
444  A(speciesNumber + 1, i) = 0;
445  }
446 }
447 
448 
449 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
450 
451 template<class CompType, class ThermoType>
453 (
454  const Foam::scalarField& phiq,
455  scalarField& Rphiq
456 )
457 {
458  bool retrieved(false);
460 
461  // If the tree is not empty
462  if (chemisTree_.size())
463  {
464  chemisTree_.binaryTreeSearch(phiq, chemisTree_.root(), phi0);
465 
466  // lastSearch keeps track of the chemPoint we obtain by the regular
467  // binary tree search
468  lastSearch_ = phi0;
469  if (phi0->inEOA(phiq))
470  {
471  retrieved = true;
472  }
473  // After a successful secondarySearch, phi0 store a pointer to the
474  // found chemPoint
475  else if (chemisTree_.secondaryBTSearch(phiq, phi0))
476  {
477  retrieved = true;
478  }
479  else if (MRURetrieve_)
480  {
481  typename SLList
482  <
484  >::iterator iter = MRUList_.begin();
485 
486  for ( ; iter != MRUList_.end(); ++iter)
487  {
488  phi0 = iter();
489  if (phi0->inEOA(phiq))
490  {
491  retrieved = true;
492  break;
493  }
494  }
495  }
496  }
497  // The tree is empty, retrieved is still false
498  else
499  {
500  // There is no chempoints that we can try to grow
501  lastSearch_ = nullptr;
502  }
503 
504  if (retrieved)
505  {
506  phi0->increaseNumRetrieve();
507  scalar elapsedTimeSteps =
508  this->chemistry_.timeSteps() - phi0->timeTag();
509 
510  // Raise a flag when the chemPoint has been used more than the allowed
511  // number of time steps
512  if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->toRemove())
513  {
514  cleaningRequired_ = true;
515  phi0->toRemove() = true;
516  }
517  lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
518  addToMRU(phi0);
519  calcNewC(phi0, phiq, Rphiq);
520  nRetrieved_++;
521  return true;
522  }
523  else
524  {
525  // This point is reached when every retrieve trials have failed
526  // or if the tree is empty
527  return false;
528  }
529 }
530 
531 
532 template<class CompType, class ThermoType>
534 (
535  const scalarField& phiq,
536  const scalarField& Rphiq,
537  const label li,
538  const scalar rho,
539  const scalar deltaT
540 )
541 {
542  label growthOrAddFlag = 1;
543  // If lastSearch_ holds a valid pointer to a chemPoint AND the growPoints_
544  // option is on, the code first tries to grow the point hold by lastSearch_
545  if (lastSearch_ && growPoints_)
546  {
547  if (grow(lastSearch_, phiq, Rphiq))
548  {
549  nGrowth_++;
550  growthOrAddFlag = 0;
551  addToMRU(lastSearch_);
552  //the structure of the tree is not modified, return false
553  return growthOrAddFlag;
554  }
555  }
556 
557  // If the code reach this point, it is either because lastSearch_ is not
558  // valid, OR because growPoints_ is not on, OR because the grow operation
559  // has failed. In the three cases, a new point is added to the tree.
560  if (chemisTree().isFull())
561  {
562  // If cleanAndBalance operation do not result in a reduction of the tree
563  // size, the last possibility is to delete completely the tree.
564  // It can be partially rebuild with the MRU list if this is used.
565  if (!cleanAndBalance())
566  {
568  if (maxMRUSize_>0)
569  {
570  // Create a copy of each chemPointISAT of the MRUList_ before
571  // they are deleted
572  typename SLList
573  <
575  >::iterator iter = MRUList_.begin();
576  for ( ; iter != MRUList_.end(); ++iter)
577  {
578  tempList.append
579  (
581  );
582  }
583  }
584  chemisTree().clear();
585 
586  // Pointers to chemPoint are not valid anymore, clear the list
587  MRUList_.clear();
588 
589  // Construct the tree without giving a reference to attach to it
590  // since the structure has been completely discarded
592  forAll(tempList, i)
593  {
594  chemisTree().insertNewLeaf
595  (
596  tempList[i]->phi(),
597  tempList[i]->Rphi(),
598  tempList[i]->A(),
599  scaleFactor(),
600  this->tolerance(),
601  scaleFactor_.size(),
602  nulPhi
603  );
604  deleteDemandDrivenData(tempList[i]);
605  }
606  }
607 
608  // The structure has been changed, it will force the binary tree to
609  // perform a new search and find the most appropriate point still stored
610  lastSearch_ = nullptr;
611  }
612 
613  // Compute the A matrix needed to store the chemPoint.
614  label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
615  scalarSquareMatrix A(ASize, Zero);
616  computeA(A, Rphiq, li, rho, deltaT);
617 
618  chemisTree().insertNewLeaf
619  (
620  phiq,
621  Rphiq,
622  A,
623  scaleFactor(),
624  this->tolerance(),
625  scaleFactor_.size(),
626  lastSearch_ // lastSearch_ may be nullptr (handled by binaryTree)
627  );
628  if (lastSearch_ != nullptr)
629  {
630  addToMRU(lastSearch_);
631  }
632  nAdd_++;
633 
634  return growthOrAddFlag;
635 }
636 
637 
638 template<class CompType, class ThermoType>
639 void
641 {
642  if (this->log())
643  {
644  nRetrievedFile_()
645  << runTime_.timeOutputValue() << " " << nRetrieved_ << endl;
646  nRetrieved_ = 0;
647 
648  nGrowthFile_()
649  << runTime_.timeOutputValue() << " " << nGrowth_ << endl;
650  nGrowth_ = 0;
651 
652  nAddFile_()
653  << runTime_.timeOutputValue() << " " << nAdd_ << endl;
654  nAdd_ = 0;
655 
656  sizeFile_()
657  << runTime_.timeOutputValue() << " " << this->size() << endl;
658  }
659 }
660 
661 
662 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
dimensionedScalar log(const dimensionedScalar &ds)
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 > &)
Class to perform the LU decomposition on a symmetric matrix.
const iterator & end()
Definition: LList.H:286
Template class for non-intrusive linked lists.
Definition: LList.H:51
T & first()
Return the first entry added.
Definition: LList.H:139
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
bool checkSolution(const scalarField &phiq, const scalarField &Rphiq)
If phiq is not in the EOA, then the mapping is computed.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Leaf of the binary tree. The chemPoint stores the composition &#39;phi&#39;, the mapping of this composition ...
List< label > & completeToSimplifiedIndex()
ISAT(const dictionary &chemistryProperties, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from dictionary.
Definition: ISAT.C:33
BasicChemistryModel< rhoReactionThermo > & chemistry
const label & timeTag()
phi
Definition: pEqn.H:104
virtual bool retrieve(const Foam::scalarField &phiq, scalarField &Rphiq)
Find the closest stored leaf of phiQ and store the result in.
Definition: ISAT.C:453
label size() const
Return the number of elements in matrix (m*n)
Definition: MatrixI.H:71
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
virtual label add(const scalarField &phiq, const scalarField &Rphiq, const label li, const scalar rho, const scalar deltaT)
Add information to the tabulation.
Definition: ISAT.C:534
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.
static const zero Zero
Definition: zero.H:97
const scalarField & phi() const
const dimensionedScalar & phi0
Magnetic flux quantum: default SI units: [Wb].
iterator begin()
Definition: LList.H:281
const scalarField & Rphi() const
An abstract class for chemistry tabulation.
virtual label nEqns() const
Number of ODE&#39;s to solve.
const Time & time() const
Return time.
Definition: IOobject.C:328
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.
autoPtr< OFstream > logFile(const word &name) const
Create and return a TDAC log file of the given name.
Implementation of the ISAT (In-situ adaptive tabulation), for chemistry calculation.
Definition: ISAT.H:57
void deleteDemandDrivenData(DataPtr &dataPtr)
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:812
void inv(scalarSquareMatrix &M) const
Set M to the inverse of this square matrix.