ISAT.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2016-2017 OpenFOAM Foundation
9 Copyright (C) 2019-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "ISAT.H"
30#include "LUscalarMatrix.H"
31#include "demandDrivenData.H"
32
33// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
34
35template<class CompType, class ThermoType>
37(
38 const dictionary& chemistryProperties,
40)
41:
42 chemistryTabulationMethod<CompType, ThermoType>
43 (
44 chemistryProperties,
46 ),
47 chemisTree_(chemistry, this->coeffsDict_),
48 scaleFactor_(chemistry.nEqns() + ((this->variableTimeStep()) ? 1 : 0), 1),
49 runTime_(chemistry.time()),
50 chPMaxLifeTime_
51 (
52 this->coeffsDict_.getOrDefault("chPMaxLifeTime", INT_MAX)
53 ),
54 maxGrowth_(this->coeffsDict_.getOrDefault("maxGrowth", INT_MAX)),
55 checkEntireTreeInterval_
56 (
57 this->coeffsDict_.getOrDefault("checkEntireTreeInterval", INT_MAX)
58 ),
59 maxDepthFactor_
60 (
61 this->coeffsDict_.getOrDefault
62 (
63 "maxDepthFactor",
64 (chemisTree_.maxNLeafs() - 1)
65 /(log(scalar(chemisTree_.maxNLeafs()))/log(2.0))
66 )
67 ),
68 minBalanceThreshold_
69 (
70 this->coeffsDict_.getOrDefault
71 (
72 "minBalanceThreshold",0.1*chemisTree_.maxNLeafs()
73 )
74 ),
75 MRURetrieve_(this->coeffsDict_.getOrDefault("MRURetrieve", false)),
76 maxMRUSize_(this->coeffsDict_.getOrDefault("maxMRUSize", 0)),
77 lastSearch_(nullptr),
78 growPoints_(this->coeffsDict_.getOrDefault("growPoints", true)),
79 nRetrieved_(0),
80 nGrowth_(0),
81 nAdd_(0),
82 cleaningRequired_(false)
83{
84 if (this->active_)
85 {
86 dictionary scaleDict(this->coeffsDict_.subDict("scaleFactor"));
87 const label Ysize = this->chemistry_.Y().size();
88 const scalar otherScaleFactor = scaleDict.get<scalar>("otherSpecies");
89 for (label i=0; i<Ysize; ++i)
90 {
91 if (!scaleDict.found(this->chemistry_.Y()[i].member()))
92 {
93 scaleFactor_[i] = otherScaleFactor;
94 }
95 else
96 {
97 scaleFactor_[i] =
98 scaleDict.get<scalar>(this->chemistry_.Y()[i].member());
99 }
100 }
101
102 scaleDict.readEntry("Temperature", scaleFactor_[Ysize]);
103 scaleDict.readEntry("Pressure", scaleFactor_[Ysize + 1]);
104
105 if (this->variableTimeStep())
106 {
107 scaleDict.readEntry("deltaT", scaleFactor_[Ysize + 2]);
108 }
109 }
110
111 if (this->variableTimeStep())
112 {
113 nAdditionalEqns_ = 3;
114 }
115 else
116 {
117 nAdditionalEqns_ = 2;
118 }
119
120 if (this->log())
121 {
122 nRetrievedFile_ = chemistry.logFile("found_isat.out");
123 nGrowthFile_ = chemistry.logFile("growth_isat.out");
124 nAddFile_ = chemistry.logFile("add_isat.out");
125 sizeFile_ = chemistry.logFile("size_isat.out");
126 }
127}
128
129
130// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
131
132template<class CompType, class ThermoType>
134{}
135
136
137// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
138
139template<class CompType, class ThermoType>
141(
143)
144{
145 if (maxMRUSize_ > 0 && MRURetrieve_)
146 {
147 auto iter = MRUList_.begin();
148
149 // First search if the chemPoint is already in the list
150 bool isInList = false;
151 for (; iter.good(); ++iter)
152 {
153 if (*iter == phi0)
154 {
155 isInList = true;
156 break;
157 }
158 }
159
160 if (isInList)
161 {
162 // If it is in the list, then move it to front
163 if (iter != MRUList_.begin())
164 {
165 MRUList_.remove(iter);
166 MRUList_.insert(phi0);
167 }
168 }
169 else
170 {
171 // chemPoint not yet in the list, iter is last
172 if (MRUList_.size() == maxMRUSize_)
173 {
174 if (iter == MRUList_.end())
175 {
176 MRUList_.remove(iter);
177 MRUList_.insert(phi0);
178 }
179 else
180 {
182 << "Error in MRUList construction"
183 << exit(FatalError);
184 }
185 }
186 else
187 {
188 MRUList_.insert(phi0);
189 }
190 }
191 }
192}
193
194
195template<class CompType, class ThermoType>
197(
198 chemPointISAT<CompType, ThermoType>* phi0,
199 const scalarField& phiq,
200 scalarField& Rphiq
201)
202{
203 label nEqns = this->chemistry_.nEqns(); // Species, T, p
204 bool mechRedActive = this->chemistry_.mechRed()->active();
205 Rphiq = phi0->Rphi();
206 scalarField dphi(phiq-phi0->phi());
207 const scalarSquareMatrix& gradientsMatrix = phi0->A();
208 List<label>& completeToSimplified(phi0->completeToSimplifiedIndex());
209
210 // Rphiq[i]=Rphi0[i]+A(i, j)dphi[j]
211 // where Aij is dRi/dphi_j
212 for (label i=0; i<nEqns-nAdditionalEqns_; ++i)
213 {
214 if (mechRedActive)
215 {
216 label si = completeToSimplified[i];
217 // The species is active
218 if (si != -1)
219 {
220 for (label j=0; j<nEqns-2; j++)
221 {
222 label sj = completeToSimplified[j];
223 if (sj != -1)
224 {
225 Rphiq[i] += gradientsMatrix(si, sj)*dphi[j];
226 }
227 }
228 Rphiq[i] +=
229 gradientsMatrix(si, phi0->nActiveSpecies())*dphi[nEqns - 2];
230 Rphiq[i] +=
231 gradientsMatrix(si, phi0->nActiveSpecies() + 1)
232 *dphi[nEqns - 1];
233
234 if (this->variableTimeStep())
235 {
236 Rphiq[i] +=
237 gradientsMatrix(si, phi0->nActiveSpecies() + 2)
238 *dphi[nEqns];
239 }
240
241 // As we use an approximation of A, Rphiq should be checked for
242 // negative values
243 Rphiq[i] = max(0.0, Rphiq[i]);
244 }
245 // The species is not active A(i, j) = I(i, j)
246 else
247 {
248 Rphiq[i] += dphi[i];
249 Rphiq[i] = max(0.0, Rphiq[i]);
250 }
251 }
252 else // Mechanism reduction is not active
253 {
254 for (label j=0; j<nEqns; ++j)
255 {
256 Rphiq[i] += gradientsMatrix(i, j)*dphi[j];
257 }
258 // As we use a first order gradient matrix, Rphiq should be checked
259 // for negative values
260 Rphiq[i] = max(0.0, Rphiq[i]);
261 }
262 }
263}
264
265
266template<class CompType, class ThermoType>
268(
269 chemPointISAT<CompType, ThermoType>* phi0,
270 const scalarField& phiq,
271 const scalarField& Rphiq
272)
273{
274 // If the pointer to the chemPoint is nullptr, the function stops
275 if (phi0 == nullptr)
276 {
277 return false;
278 }
279
280 // Raise a flag when the chemPoint used has been grown more than the
281 // allowed number of time
282 if (phi0->nGrowth() > maxGrowth_)
283 {
284 cleaningRequired_ = true;
285 phi0->toRemove() = true;
286 return false;
287 }
288
289 // If the solution RphiQ is still within the tolerance we try to grow it
290 // in some cases this might result in a failure and the grow function of
291 // the chemPoint returns false
292 if (phi0->checkSolution(phiq, Rphiq))
293 {
294 return phi0->grow(phiq);
295 }
296
297 // The actual solution and the approximation given by ISAT are too different
298 return false;
299}
300
301
302template<class CompType, class ThermoType>
303bool
305{
306 bool treeModified(false);
307
308 // Check all chemPoints to see if we need to delete some of the chemPoints
309 // according to the elapsed time and number of growths
310 chemPointISAT<CompType, ThermoType>* x = chemisTree_.treeMin();
311 while (x != nullptr)
312 {
313 chemPointISAT<CompType, ThermoType>* xtmp =
314 chemisTree_.treeSuccessor(x);
315
316 scalar elapsedTimeSteps = this->chemistry_.timeSteps() - x->timeTag();
317
318 if ((elapsedTimeSteps > chPMaxLifeTime_) || (x->nGrowth() > maxGrowth_))
319 {
320 chemisTree_.deleteLeaf(x);
321 treeModified = true;
322 }
323 x = xtmp;
324 }
325 // Check if the tree should be balanced according to criterion:
326 // - the depth of the tree bigger than a*log2(size), log2(size) being the
327 // ideal depth (e.g. 4 leafs can be stored in a tree of depth 2)
328 if
329 (
330 chemisTree_.size() > minBalanceThreshold_
331 && chemisTree_.depth() >
332 maxDepthFactor_*log(scalar(chemisTree_.size()))/log(2.0)
333 )
334 {
335 chemisTree_.balance();
336 MRUList_.clear();
337 treeModified = true;
338 }
339
340 // Return a bool to specify if the tree structure has been modified and is
341 // now below the user specified limit (true if not full)
342 return (treeModified && !chemisTree_.isFull());
343}
344
345
346template<class CompType, class ThermoType>
348(
350 const scalarField& Rphiq,
351 const scalar rhoi,
352 const scalar dt
353)
354{
355 bool mechRedActive = this->chemistry_.mechRed()->active();
356 label speciesNumber = this->chemistry_.nSpecie();
357 scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
358 for (label i=0; i<speciesNumber; ++i)
359 {
360 label s2c = i;
361 if (mechRedActive)
362 {
363 s2c = this->chemistry_.simplifiedToCompleteIndex()[i];
364 }
365 Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
366 }
367 Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
368 Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
369 if (this->variableTimeStep())
370 {
371 Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
372 }
373
374 // Aaa is computed implicitly,
375 // A is given by A = C(psi0, t0+dt), where C is obtained through solving
376 // d/dt C(psi0,t) = J(psi(t))C(psi0,t)
377 // If we solve it implicitly:
378 // (C(psi0, t0+dt) - C(psi0,t0))/dt = J(psi(t0+dt))C(psi0,t0+dt)
379 // The Jacobian is thus computed according to the mapping
380 // C(psi0,t0+dt)*(I-dt*J(psi(t0+dt))) = C(psi0, t0)
381 // A = C(psi0,t0)/(I-dt*J(psi(t0+dt)))
382 // where C(psi0,t0) = I
383
384 this->chemistry_.jacobian(runTime_.value(), Rcq, A);
385
386 // The jacobian is computed according to the molar concentration
387 // the following conversion allows the code to use A with mass fraction
388 for (label i=0; i<speciesNumber; ++i)
389 {
390 label si = i;
391
392 if (mechRedActive)
393 {
394 si = this->chemistry_.simplifiedToCompleteIndex()[i];
395 }
396
397 for (label j=0; j<speciesNumber; ++j)
398 {
399 label sj = j;
400 if (mechRedActive)
401 {
402 sj = this->chemistry_.simplifiedToCompleteIndex()[j];
403 }
404 A(i, j) *=
405 -dt*this->chemistry_.specieThermo()[si].W()
406 /this->chemistry_.specieThermo()[sj].W();
407 }
408
409 A(i, i) += 1;
410 // Columns for pressure and temperature
411 A(i, speciesNumber) *=
412 -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
413 A(i, speciesNumber+1) *=
414 -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
415 }
416
417 // For temperature and pressure, only unity on the diagonal
418 A(speciesNumber, speciesNumber) = 1;
419 A(speciesNumber + 1, speciesNumber + 1) = 1;
420 if (this->variableTimeStep())
421 {
422 A[speciesNumber + 2][speciesNumber + 2] = 1;
423 }
424
425 // Inverse of (I-dt*J(psi(t0+dt)))
426 LUscalarMatrix LUA(A);
427 LUA.inv(A);
428}
429
430
431// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
432
433template<class CompType, class ThermoType>
435(
436 const Foam::scalarField& phiq,
437 scalarField& Rphiq
438)
439{
440 bool retrieved(false);
442
443 // If the tree is not empty
444 if (chemisTree_.size())
445 {
446 chemisTree_.binaryTreeSearch(phiq, chemisTree_.root(), phi0);
447
448 // lastSearch keeps track of the chemPoint we obtain by the regular
449 // binary tree search
450 lastSearch_ = phi0;
451 if (phi0->inEOA(phiq))
452 {
453 retrieved = true;
454 }
455 // After a successful secondarySearch, phi0 store a pointer to the
456 // found chemPoint
457 else if (chemisTree_.secondaryBTSearch(phiq, phi0))
458 {
459 retrieved = true;
460 }
461 else if (MRURetrieve_)
462 {
463 forAllConstIters(MRUList_, iter)
464 {
465 phi0 = iter();
466 if (phi0->inEOA(phiq))
467 {
468 retrieved = true;
469 break;
470 }
471 }
472 }
473 }
474 // The tree is empty, retrieved is still false
475 else
476 {
477 // There is no chempoints that we can try to grow
478 lastSearch_ = nullptr;
479 }
480
481 if (retrieved)
482 {
483 phi0->increaseNumRetrieve();
484 scalar elapsedTimeSteps =
485 this->chemistry_.timeSteps() - phi0->timeTag();
486
487 // Raise a flag when the chemPoint has been used more than the allowed
488 // number of time steps
489 if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->toRemove())
490 {
491 cleaningRequired_ = true;
492 phi0->toRemove() = true;
493 }
494 lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
495 addToMRU(phi0);
496 calcNewC(phi0,phiq, Rphiq);
497 ++nRetrieved_;
498 return true;
499 }
500
501
502 // This point is reached when every retrieve trials have failed
503 // or if the tree is empty
504 return false;
505}
506
507
508template<class CompType, class ThermoType>
510(
511 const scalarField& phiq,
512 const scalarField& Rphiq,
513 const scalar rho,
514 const scalar deltaT
515)
516{
517 label growthOrAddFlag = 1;
518
519 // If lastSearch_ holds a valid pointer to a chemPoint AND the growPoints_
520 // option is on, the code first tries to grow the point hold by lastSearch_
521 if (lastSearch_ && growPoints_)
522 {
523 if (grow(lastSearch_,phiq, Rphiq))
524 {
525 ++nGrowth_;
526 growthOrAddFlag = 0;
527
528 // The structure of the tree is not modified, return false
529 return growthOrAddFlag;
530 }
531 }
532
533 // If the code reach this point, it is either because lastSearch_ is not
534 // valid, OR because growPoints_ is not on, OR because the grow operation
535 // has failed. In the three cases, a new point is added to the tree.
536 if (chemisTree().isFull())
537 {
538 // If cleanAndBalance operation do not result in a reduction of the tree
539 // size, the last possibility is to delete completely the tree.
540 // It can be partially rebuilt with the MRU list if this is used.
541 if (!cleanAndBalance())
542 {
544 if (maxMRUSize_>0)
545 {
546 // Create a copy of each chemPointISAT of the MRUList_ before
547 // they are deleted
548 forAllConstIters(MRUList_, iter)
549 {
550 tempList.append
551 (
553 );
554 }
555 }
556 chemisTree().clear();
557
558 // Pointers to chemPoint are not valid anymore, clear the list
559 MRUList_.clear();
560
561 // Construct the tree without giving a reference to attach to it
562 // since the structure has been completely discarded
563 chemPointISAT<CompType, ThermoType>* nulPhi = nullptr;
564 for (auto& t : tempList)
565 {
566 chemisTree().insertNewLeaf
567 (
568 t->phi(),
569 t->Rphi(),
570 t->A(),
571 scaleFactor(),
572 this->tolerance(),
573 scaleFactor_.size(),
574 nulPhi
575 );
577 }
578 }
579
580 // The structure has been changed, it will force the binary tree to
581 // perform a new search and find the most appropriate point still stored
582 lastSearch_ = nullptr;
583 }
584
585 // Compute the A matrix needed to store the chemPoint.
586 label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
587 scalarSquareMatrix A(ASize, Zero);
588 computeA(A, Rphiq, rho, deltaT);
589
590 chemisTree().insertNewLeaf
591 (
592 phiq,
593 Rphiq,
594 A,
595 scaleFactor(),
596 this->tolerance(),
597 scaleFactor_.size(),
598 lastSearch_ // lastSearch_ may be nullptr (handled by binaryTree)
599 );
600
601 ++nAdd_;
602
603 return growthOrAddFlag;
604}
605
606
607template<class CompType, class ThermoType>
608void
610{
611 if (this->log())
612 {
613 nRetrievedFile_()
614 << runTime_.timeOutputValue() << " " << nRetrieved_ << endl;
615 nRetrieved_ = 0;
616
617 nGrowthFile_()
618 << runTime_.timeOutputValue() << " " << nGrowth_ << endl;
619 nGrowth_ = 0;
620
621 nAddFile_()
622 << runTime_.timeOutputValue() << " " << nAdd_ << endl;
623 nAdd_ = 0;
624
625 sizeFile_()
626 << runTime_.timeOutputValue() << " " << this->size() << endl;
627 }
628}
629
630
631// ************************************************************************* //
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
Extends StandardChemistryModel by adding the TDAC method.
PtrList< volScalarField > & Y()
Leaf of the binary tree. The chemPoint stores the composition 'phi', the mapping of this composition ...
An abstract class for chemistry tabulation.
Switch active_
Is tabulation active?
TDACChemistryModel< CompType, ThermoType > & chemistry_
Implementation of the ISAT (In-situ adaptive tabulation), for chemistry calculation.
Definition: ISAT.H:62
virtual bool retrieve(const Foam::scalarField &phiq, scalarField &Rphiq)
Definition: ISAT.C:435
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
Definition: dictionaryI.H:87
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
Sums a given list of (at least two or more) fields and outputs the result into a new field,...
Definition: add.H:161
BasicChemistryModel< psiReactionThermo > & chemistry
Template functions to aid in the implementation of demand driven data.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const dimensionedScalar phi0
Magnetic flux quantum: default SI units: [Wb].
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar log(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
SquareMatrix< scalar > scalarSquareMatrix
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
void deleteDemandDrivenData(DataPtr &dataPtr)
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278