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-2022 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 #include "scalarField.H"
126 #include "LUscalarMatrix.H"
127 #include "Switch.H"
128 
129 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
130 
131 namespace Foam
132 {
133 
134 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
135 
136 class binaryNode;
138 namespace chemistryTabulationMethods
139 {
140  class ISAT;
141 }
142 
143 
144 /*---------------------------------------------------------------------------*\
145  Class chemPointISAT Declaration
146 \*---------------------------------------------------------------------------*/
148 class chemPointISAT
149 {
150  // Private Data
151 
152  //- Reference to the ISAT table
154 
155  //- Vector storing the composition, temperature and pressure
156  // and deltaT if a variable time step is set on
157  scalarField phi_;
158 
159  //- Vector storing the mapping of the composition phi
160  scalarField Rphi_;
161 
162  //- LT the transpose of the L matrix describing the Ellipsoid Of
163  // Accuracy use List of Lists to be able to change size if DAC is used
164  scalarSquareMatrix LT_;
165 
166  //- A the mapping gradient matrix
168 
169  //- Vector storing the scale factor
170  scalarField scaleFactor_;
171 
172  //- Reference to the node in the binary tree holding this chemPoint
173  binaryNode* node_;
174 
175  //- The size of the composition space (size of the vector phi)
176  label completeSpaceSize_;
177 
178  //- Number of times the element has been grown
179  label nGrowth_;
180 
181  //- Tolerance for the Ellipsoid of accuracy
182  static scalar tolerance_;
183 
184  //- Number of active species stored in the chemPoint
185  label nActive_;
186 
187  //- Vectors that store the index conversion between the simplified
188  // and the complete chemical mechanism
189  List<label> simplifiedToCompleteIndex_;
190 
191  label timeTag_;
192  label lastTimeUsed_;
193 
194  bool toRemove_;
195 
196  label maxNumNewDim_;
197 
198  Switch printProportion_;
199 
200  //- Variable to store the number of retrieves the chemPoint
201  // will generate at each time step
202  label numRetrieve_;
203 
204  //- Variable to store the number of time steps the chempoint is allowed
205  // to still live according to the maxChPLifeTime_ parameter
206  label nLifeTime_;
207 
208  List<label> completeToSimplifiedIndex_;
209 
210  label idT_;
211  label idp_;
212  label iddeltaT_;
213 
214  //- QR decomposition of a matrix
215  // Input : nCols cols number
216  // R the matrix to decompose
217  // QT an empty matrix that stores the transpose of the Q matrix
218  // R is returned in the given R matrix
219  // which is used to store the ellipsoid of accuracy
220  void qrDecompose
221  (
222  const label nCols,
224  );
225 
226  //- QR update of the matrix A
227  void qrUpdate
228  (
230  const label n,
231  const scalarField& u,
232  const scalarField& v
233  );
234 
235  void rotate
236  (
238  const label i,
239  const scalar a,
240  const scalar b,
241  label n
242  );
243 
244 
245 public:
246 
247  // Constructors
248 
249  //- Construct from components
251  (
253  const scalarField& phi,
254  const scalarField& Rphi,
255  const scalarSquareMatrix& A,
256  const scalarField& scaleFactor,
257  const scalar tolerance,
258  const label completeSpaceSize,
259  const label nActive,
260  const dictionary& coeffsDict,
261  binaryNode* node = nullptr
262  );
263 
264  //- Construct from another chemPoint and reference to a binary node
266  (
267  const chemPointISAT& p,
268  binaryNode* node
269  );
270 
271  //- Construct from another chemPoint
273 
274 
275  // Member Functions
276 
277  //- Access to the ISAT table
278  inline const chemistryTabulationMethods::ISAT& table() const;
279 
280  inline label nGrowth() const;
281 
282  inline label& completeSpaceSize();
283 
284  inline const label& completeSpaceSize() const;
285 
286  inline const scalarField& phi() const;
287 
288  inline const scalarField& Rphi() const;
289 
290  inline const scalarField& scaleFactor() const;
291 
292  inline const scalar& tolerance() const;
293 
294  inline static void changeTolerance(scalar newTol);
295 
296  inline binaryNode*& node();
297 
298  inline const scalarSquareMatrix& A() const;
299 
300  inline scalarSquareMatrix& A();
301 
302  inline const scalarSquareMatrix& LT() const;
303 
304  inline scalarSquareMatrix& LT();
305 
306  inline label nActive() const;
307 
308  inline const List<label>& completeToSimplifiedIndex() const;
309 
310  inline const List<label>& simplifiedToCompleteIndex() const;
311 
312  //- Increases the number of retrieves the chempoint has generated
313  inline void increaseNumRetrieve();
314 
315  //- Resets the number of retrieves at each time step
316  inline void resetNumRetrieve();
317 
318  //- Increases the "counter" of the chP life
319  inline void increaseNLifeTime();
320 
321  inline label simplifiedToCompleteIndex(const label i);
322 
323  inline const label& timeTag();
324 
325  inline label& lastTimeUsed();
326 
327  inline bool& toRemove();
328 
329  inline label& maxNumNewDim();
330 
331  inline const label& numRetrieve();
332 
333  inline const label& nLifeTime();
334 
335  // ISAT functions
336 
337  //- To RETRIEVE the mapping from the stored chemPoint phi, the query
338  // point phiq has to be in the EOA of phi.
339  // To test if phiq is in the ellipsoid:
340  // ||L^T.dphi|| <= 1
341  bool inEOA(const scalarField& phiq);
342 
343  //- More details about the minimum-volume ellipsoid covering an
344  // ellipsoid E and a point p are found in [1].
345  // Here is the main steps to obtain the
346  // modified matrix L' describing the new ellipsoid.
347  // 1) calculate the point p' in the transformed space :
348  // p' = L^T.(p-phi)
349  // 2) compute the rank-one decomposition:
350  // G = I + gamma.p'.p'^T,
351  // with gamma = (1/|p'|-1)*1/|p'|^2
352  // 3) compute L':
353  // L'L'^T = (L.G)(L.G)^T,
354  // L'^T is then obtained by QR decomposition of (L.G)^T = G^T.L^T
355  // [1] Stephen B. Pope, "Algorithms for ellipsoids", FDA 08-01,
356  // Cornell University, 2008
357  bool grow(const scalarField& phiq);
358 
359  //- If phiq is not in the EOA, then the mapping is computed.
360  // But as the EOA is a conservative approximation of the region of
361  // accuracy surrounding the point phi, we could expand it by
362  // comparing the computed results with the one obtained by linear
363  // interpolation. The error eps is calculated:
364  // eps = ||B.(dR - dRl)||,
365  // with dR = Rphiq - Rphi, dRl = A.dphi and B the diagonal scale
366  // factor matrix.
367  // If eps <= tolerance, the EOA is too conservative and a GROW is
368  // performed,
369  // otherwise, the newly computed mapping is associated to the
370  // initial composition and added to the tree.
371  bool checkSolution
372  (
373  const scalarField& phiq,
374  const scalarField& Rphiq
375  );
376 };
377 
378 
379 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
380 
381 } // End namespace Foam
382 
383 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
384 
385 #include "chemPointISATI.H"
386 
387 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
388 
389 #endif
390 
391 // ************************************************************************* //
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
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 ...
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
phi
Definition: correctPhi.H:3
Node of the binary tree.
Definition: binaryNode.H:48
#define R(A, B, C, D, E, F, K, M)
label n
Implementation of the ISAT (In-situ adaptive tabulation), for chemistry calculation.
Definition: ISAT.H:60
volScalarField & p
Namespace for OpenFOAM.