All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
chemPointISAT.H
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-2021 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 Class
25  Foam::chemPointISAT
26 
27 Description
28  Leaf of the binary tree.
29  The chemPoint stores the composition 'phi', the mapping of this
30  composition Rphi, the mapping gradient matrix A and the matrix describing
31  the Ellipsoid Of Accuracy (EOA).
32 
33  1) When the chemPoint is created the region of accuracy is approximated by
34  an ellipsoid E centered in 'phi' (obtained with the constant):
35  E = {x| ||L^T.(x-phi)|| <= 1},
36  with x a point in the composition space and L^T the transpose of an
37  upper triangular matrix describing the EOA (see below: "Computation of
38  L" ).
39 
40  2) To RETRIEVE the mapping from the chemPoint phi, the query point phiq has
41  to be in the EOA of phi. It follows that, dphi=phiq-phi and to test if
42  phiq is in the ellipsoid there are two methods. First, compare
43  r=||dphi|| with rmin and rmax. If r < rmin, phiq is in the EOA. If r >
44  rmax, phiq is out of the EOA. This operations is O(completeSpaceSize)
45  and is performed first. If rmin < r < rmax, then the second method is
46  used:
47  ||L^T.dphi|| <= 1 to be in the EOA.
48  If phiq is in the EOA, Rphiq is obtained by linear interpolation:
49  Rphiq= Rphi + A.dphi.
50 
51  3) If phiq is not in the EOA, then the mapping is computed. But as the EOA
52  is a conservative approximation of the region of accuracy surrounding
53  the point phi, we could expand it by comparing the computed results with
54  the one obtained by linear interpolation. The error epsGrow is
55  calculated:
56  epsGrow = ||B.(dR - dRl)||,
57  with dR = Rphiq - Rphi, dRl = A.dphi and B the diagonal scale factor
58  matrix. If epsGrow <= tolerance, the EOA is too conservative and a GROW
59  is performed otherwise, the newly computed mapping is associated to the
60  initial composition and added to the tree.
61 
62  4) To GROW the EOA, we expand it to include the previous EOA and the query
63  point phiq. The rank-one matrix method is used. The EOA is transformed
64  to a hypersphere centered at the origin. Then it is expanded to include
65  the transformed point phiq' on its boundary. Then the inverse
66  transformation give the modified matrix L' (see below: "Grow the EOA").
67 
68 
69  Computation of L :
70  In [1], the EOA of the constant approximation is given by
71  E = {x| ||B.A/tolerance.(x-phi)|| <= 1},
72 
73  with B a scale factor diagonal matrix, A the mapping gradient matrix and
74  tolerance the absolute tolerance. If we take the QR decomposition of
75  (B.A)/tolerance= Q.R, with Q an orthogonal matrix and R an upper
76  triangular matrix such that the EOA is described by
77  (phiq-phi0)^T.R^T.R.(phiq-phi0) <= 1
78  L^T = R, both Cholesky decomposition of A^T.B^T.B.A/tolerance^2
79  This representation of the ellipsoid is used in [2] and in order to avoid
80  large value of semi-axe length in certain direction, a Singular Value
81  Decomposition (SVD) is performed on the L matrix:
82  L = UDV^T,
83  with the orthogonal matrix U giving the directions of the principal axes
84  and 1/di the inverse of the element of the diagonal matrix D giving the
85  length of the principal semi-axes. To avoid very large value of those
86  length, di' = max(di, 1/(alphaEOA*sqrt(tolerance))), with alphaEOA = 0.1
87  (see [2]) di' = max(di, 1/2), see [1]. The latter will be used in this
88  implementation. And L' = UD'V^T, with D' the diagonal matrix with the
89  modified di'.
90 
91  Grow the EOA :
92  More details about the minimum-volume ellipsoid covering an ellipsoid E
93  and a point p are found in [3]. Here is the main steps to obtain the
94  modified matrix L' describing the new ellipsoid.
95 
96  1) calculate the point p' in the transformed space :
97  p' = L^T.(p-phi)
98  2) compute the rank-one decomposition:
99  G = I + gamma.p'.p'^T,
100  with gamma = (1/|p'|-1)*1/|p'|^2
101  3) compute L':
102  L' = L.G.
103 
104  References:
105  \verbatim
106  [1] Pope, S. B. (1997).
107  Computationally efficient implementation of combustion chemistry using
108  in situ adaptive tabulation.
109  Combustion Theory and Modelling, 1, 41-63.
110 
111  [2] Lu, L., & Pope, S. B. (2009).
112  An improved algorithm for in situ adaptive tabulation.
113  Journal of Computational Physics, 228(2), 361-386.
114 
115  [3] Pope, S. B. (2008).
116  Algorithms for ellipsoids.
117  Cornell University Report No. FDA, 08-01.
118  \endverbatim
119 
120 \*---------------------------------------------------------------------------*/
121 
122 #ifndef chemPointISAT_H
123 #define chemPointISAT_H
124 
125 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
126 
127 namespace Foam
128 {
129 
130 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
131 
132 template<class ThermoType>
133 class binaryNode;
134 
135 template<class ThermoType>
136 class TDACChemistryModel;
137 
138 
139 /*---------------------------------------------------------------------------*\
140  Class chemPointISAT Declaration
141 \*---------------------------------------------------------------------------*/
142 
143 template<class ThermoType>
144 class chemPointISAT
145 {
146  // Private Data
147 
148  //- Pointer to the chemistryModel object
149  TDACChemistryModel<ThermoType>& chemistry_;
150 
151  //- Vector storing the composition, temperature and pressure
152  // and deltaT if a variable time step is set on
153  scalarField phi_;
154 
155  //- Vector storing the mapping of the composition phi
156  scalarField Rphi_;
157 
158  //- LT the transpose of the L matrix describing the Ellipsoid Of
159  // Accuracy use List of Lists to be able to change size if DAC is used
160  scalarSquareMatrix LT_;
161 
162  //- A the mapping gradient matrix
164 
165  //- Vector storing the scale factor
166  scalarField scaleFactor_;
167 
168  //- Reference to the node in the binary tree holding this chemPoint
169  binaryNode<ThermoType>* node_;
170 
171  //- The size of the composition space (size of the vector phi)
172  label completeSpaceSize_;
173 
174  //- Number of times the element has been grown
175  label nGrowth_;
176 
177  //- Tolerance for the Ellipsoid of accuracy
178  static scalar tolerance_;
179 
180  //- Number of active species stored in the chemPoint
181  label nActiveSpecies_;
182 
183  //- Vectors that store the index conversion between the simplified
184  // and the complete chemical mechanism
185  List<label> simplifiedToCompleteIndex_;
186 
187  label timeTag_;
188  label lastTimeUsed_;
189 
190  bool toRemove_;
191 
192  label maxNumNewDim_;
193 
194  Switch printProportion_;
195 
196  //- Variable to store the number of retrieves the chemPoint
197  // will generate at each time step
198  label numRetrieve_;
199 
200  //- Variable to store the number of time steps the chempoint is allowed
201  // to still live according to the maxChPLifeTime_ parameter
202  label nLifeTime_;
203 
204  List<label> completeToSimplifiedIndex_;
205 
206  label idT_;
207  label idp_;
208  label iddeltaT_;
209 
210  //- QR decomposition of a matrix
211  // Input : nCols cols number
212  // R the matrix to decompose
213  // QT an empty matrix that stores the transpose of the Q matrix
214  // R is returned in the given R matrix
215  // which is used to store the ellipsoid of accuracy
216  void qrDecompose
217  (
218  const label nCols,
220  );
221 
222  //- QR update of the matrix A
223  void qrUpdate
224  (
226  const label n,
227  const scalarField& u,
228  const scalarField& v
229  );
230 
231  void rotate
232  (
234  const label i,
235  const scalar a,
236  const scalar b,
237  label n
238  );
239 
240 
241 public:
242 
243  // Constructors
244 
245  //- Construct from components
247  (
249  const scalarField& phi,
250  const scalarField& Rphi,
251  const scalarSquareMatrix& A,
252  const scalarField& scaleFactor,
253  const scalar& tolerance,
254  const label& completeSpaceSize,
255  const dictionary& coeffsDict,
256  binaryNode<ThermoType>* node = nullptr
257  );
258 
259  //- Construct from another chemPoint and reference to a binary node
261  (
264  );
265 
266  //- Construct from another chemPoint
268 
269 
270  // Member Functions
271 
272  //- Access to the TDACChemistryModel
274  {
275  return chemistry_;
276  }
278  inline label nGrowth()
279  {
280  return nGrowth_;
281  }
283  inline label& completeSpaceSize()
284  {
285  return completeSpaceSize_;
286  }
288  inline const scalarField& phi() const
289  {
290  return phi_;
291  }
293  inline const scalarField& Rphi() const
294  {
295  return Rphi_;
296  }
298  inline const scalarField& scaleFactor()
299  {
300  return scaleFactor_;
301  }
303  inline const scalar& tolerance()
304  {
305  return tolerance_;
306  }
308  static void changeTolerance(scalar newTol)
309  {
310  tolerance_ = newTol;
311  }
313  inline binaryNode<ThermoType>*& node()
314  {
315  return node_;
316  }
318  inline const scalarSquareMatrix& A() const
319  {
320  return A_;
321  }
323  inline scalarSquareMatrix& A()
324  {
325  return A_;
326  }
328  inline const scalarSquareMatrix& LT() const
329  {
330  return LT_;
331  }
333  inline scalarSquareMatrix& LT()
334  {
335  return LT_;
336  }
338  inline label nActiveSpecies()
339  {
340  return nActiveSpecies_;
341  }
344  {
345  return completeToSimplifiedIndex_;
346  }
349  {
350  return simplifiedToCompleteIndex_;
351  }
352 
353  //- Increases the number of retrieves the chempoint has generated
354  void increaseNumRetrieve();
355 
356  //- Resets the number of retrieves at each time step
357  void resetNumRetrieve();
358 
359  //- Increases the "counter" of the chP life
360  void increaseNLifeTime();
361 
364  inline const label& timeTag()
365  {
366  return timeTag_;
367  }
369  inline label& lastTimeUsed()
370  {
371  return lastTimeUsed_;
372  }
374  inline bool& toRemove()
375  {
376  return toRemove_;
377  }
379  inline label& maxNumNewDim()
380  {
381  return maxNumNewDim_;
382  }
384  inline const label& numRetrieve()
385  {
386  return numRetrieve_;
387  }
389  inline const label& nLifeTime()
390  {
391  return nLifeTime_;
392  }
393 
394  // ISAT functions
395 
396  //- To RETRIEVE the mapping from the stored chemPoint phi, the query
397  // point phiq has to be in the EOA of phi.
398  // To test if phiq is in the ellipsoid:
399  // ||L^T.dphi|| <= 1
400  bool inEOA(const scalarField& phiq);
401 
402  //- More details about the minimum-volume ellipsoid covering an
403  // ellipsoid E and a point p are found in [1].
404  // Here is the main steps to obtain the
405  // modified matrix L' describing the new ellipsoid.
406  // 1) calculate the point p' in the transformed space :
407  // p' = L^T.(p-phi)
408  // 2) compute the rank-one decomposition:
409  // G = I + gamma.p'.p'^T,
410  // with gamma = (1/|p'|-1)*1/|p'|^2
411  // 3) compute L':
412  // L'L'^T = (L.G)(L.G)^T,
413  // L'^T is then obtained by QR decomposition of (L.G)^T = G^T.L^T
414  // [1] Stephen B. Pope, "Algorithms for ellipsoids", FDA 08-01,
415  // Cornell University, 2008
416  bool grow(const scalarField& phiq);
417 
418  //- If phiq is not in the EOA, then the mapping is computed.
419  // But as the EOA is a conservative approximation of the region of
420  // accuracy surrounding the point phi, we could expand it by
421  // comparing the computed results with the one obtained by linear
422  // interpolation. The error eps is calculated:
423  // eps = ||B.(dR - dRl)||,
424  // with dR = Rphiq - Rphi, dRl = A.dphi and B the diagonal scale
425  // factor matrix.
426  // If eps <= tolerance, the EOA is too conservative and a GROW is
427  // performed,
428  // otherwise, the newly computed mapping is associated to the
429  // initial composition and added to the tree.
430  bool checkSolution
431  (
432  const scalarField& phiq,
433  const scalarField& Rphiq
434  );
435 };
436 
437 
438 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
439 
440 } // End namespace Foam
441 
442 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
443 
444 #ifdef NoRepository
445  #include "chemPointISAT.C"
446 #endif
447 
448 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
449 
450 #endif
451 
452 // ************************************************************************* //
Extends standardChemistryModel by adding the TDAC method.
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
List< label > & simplifiedToCompleteIndex()
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
bool checkSolution(const scalarField &phiq, const scalarField &Rphiq)
If phiq is not in the EOA, then the mapping is computed.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Leaf of the binary tree. The chemPoint stores the composition &#39;phi&#39;, the mapping of this composition ...
List< label > & completeToSimplifiedIndex()
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none/any.
Definition: Switch.H:60
const label & timeTag()
chemPointISAT(TDACChemistryModel< ThermoType > &chemistry, const scalarField &phi, const scalarField &Rphi, const scalarSquareMatrix &A, const scalarField &scaleFactor, const scalar &tolerance, const label &completeSpaceSize, const dictionary &coeffsDict, binaryNode< ThermoType > *node=nullptr)
Construct from components.
binaryNode< ThermoType > *& node()
const label & nLifeTime()
const scalarField & scaleFactor()
TDACChemistryModel< ThermoType > & chemistry()
Access to the TDACChemistryModel.
const label & numRetrieve()
bool grow(const scalarField &phiq)
More details about the minimum-volume ellipsoid covering an.
Node of the binary tree.
Definition: binaryNode.H:49
const scalarField & phi() const
const scalarField & Rphi() const
void resetNumRetrieve()
Resets the number of retrieves at each time step.
#define R(A, B, C, D, E, F, K, M)
void increaseNLifeTime()
Increases the "counter" of the chP life.
label & completeSpaceSize()
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.
label n
volScalarField & p
static void changeTolerance(scalar newTol)
Namespace for OpenFOAM.