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-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 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 upper
37  triangular matrix describing the EOA (see below: "Computation of L" ).
38 
39  2)To RETRIEVE the mapping from the chemPoint phi, the query point phiq has to
40  be in the EOA of phi. It follows that, dphi=phiq-phi and to test if phiq
41  is in the ellipsoid there are two methods. First, compare r=||dphi|| with
42  rmin and rmax. If r < rmin, phiq is in the EOA. If r > rmax, phiq is out of
43  the EOA. This operations is O(completeSpaceSize) and is performed first.
44  If rmin < r < rmax, then the second method is used:
45  ||L^T.dphi|| <= 1 to be in the EOA.
46 
47  If phiq is in the EOA, Rphiq is obtained by linear interpolation:
48  Rphiq= Rphi + A.dphi.
49 
50  3)If phiq is not in the EOA, then the mapping is computed. But as the EOA
51  is a conservative approximation of the region of accuracy surrounding the
52  point phi, we could expand it by comparing the computed results with the
53  one obtained by linear interpolation. The error epsGrow is calculated:
54  epsGrow = ||B.(dR - dRl)||,
55  with dR = Rphiq - Rphi, dRl = A.dphi and B the diagonal scale factor
56  matrix.
57  If epsGrow <= tolerance, the EOA is too conservative and a GROW is perforned
58  otherwise, the newly computed mapping is associated to the initial
59  composition and added to the tree.
60 
61  4)To GROW the EOA, we expand it to include the previous EOA and the query
62  point phiq. The rank-one matrix method is used. The EOA is transformed
63  to a hypersphere centered at the origin. Then it is expanded to include
64  the transformed point phiq' on its boundary. Then the inverse transformation
65  give the modified matrix L' (see below: "Grow the EOA").
66 
67 
68  Computation of L :
69  In [1], the EOA of the constant approximation is given by
70  E = {x| ||B.A/tolerance.(x-phi)|| <= 1},
71  with B a scale factor diagonal matrix, A the mapping gradient matrix and
72  tolerance the absolute tolerance. If we take the QR decomposition of
73  (B.A)/tolerance= Q.R, with Q an orthogonal matrix and R an upper triangular
74  matrix such that the EOA is described by
75  (phiq-phi0)^T.R^T.R.(phiq-phi0) <= 1
76  L^T = R, both Cholesky decomposition of A^T.B^T.B.A/tolerance^2
77  This representation of the ellipsoid is used in [2] and in order to avoid
78  large value of semi-axe length in certain direction, a Singular Value
79  Decomposition (SVD) is performed on the L matrix:
80  L = UDV^T,
81  with the orthogonal matrix U giving the directions of the principal axes
82  and 1/di the inverse of the element of the diagonal matrix D giving the
83  length of the principal semi-axes. To avoid very large value of those
84  length,
85  di' = max(di, 1/(alphaEOA*sqrt(tolerance))), with alphaEOA = 0.1 (see [2])
86  di' = max(di, 1/2), see [1]. The latter will be used in this implementation.
87  And L' = UD'V^T, with D' the diagonal matrix with the modified di'.
88 
89  Grow the EOA :
90  More details about the minimum-volume ellipsoid covering an ellispoid E and
91  a point p are found in [3]. Here is the main steps to obtain the modified
92  matrix L' describind the new ellipsoid.
93  1) calculate the point p' in the transformed space :
94  p' = L^T.(p-phi)
95  2) compute the rank-one decomposition:
96  G = I + gamma.p'.p'^T,
97  with gamma = (1/|p'|-1)*1/|p'|^2
98  3) compute L':
99  L' = L.G.
100 
101  References:
102  \verbatim
103  [1] Pope, S. B. (1997).
104  Computationally efficient implementation of combustion chemistry using
105  in situ adaptive tabulation.
106  Combustion Theory and Modelling, 1, 41-63.
107 
108  [2] Lu, L., & Pope, S. B. (2009).
109  An improved algorithm for in situ adaptive tabulation.
110  Journal of Computational Physics, 228(2), 361-386.
111 
112  [3] Pope, S. B. (2008).
113  Algorithms for ellipsoids.
114  Cornell University Report No. FDA, 08-01.
115  \endverbatim
116 
117 \*---------------------------------------------------------------------------*/
118 
119 #ifndef chemPointISAT_H
120 #define chemPointISAT_H
121 
122 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
123 
124 namespace Foam
125 {
126 
127 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
128 
129 template<class CompType, class ThermoType>
130 class binaryNode;
131 
132 template<class CompType, class ThermoType>
133 class TDACChemistryModel;
134 
135 
136 /*---------------------------------------------------------------------------*\
137  Class chemPointISAT Declaration
138 \*---------------------------------------------------------------------------*/
139 
140 template<class CompType, class ThermoType>
141 class chemPointISAT
142 {
143  // Private data
144 
145  //- Pointer to the chemistryModel object
147 
148  //- Vector storing the composition, temperature and pressure
149  // and deltaT if a variable time step is set on
150  scalarField phi_;
151 
152  //- Vector storing the mapping of the composition phi
153  scalarField Rphi_;
154 
155  //- LT the transpose of the L matrix describing the Ellipsoid Of
156  // Accuracy use List of Lists to be able to change size if DAC is used
157  scalarSquareMatrix LT_;
158 
159  //- A the mapping gradient matrix
161 
162  //- Vector storing the scale factor
163  scalarField scaleFactor_;
164 
165  //- Reference to the node in the binary tree holding this chemPoint
167 
168  //- The size of the composition space (size of the vector phi)
169  label completeSpaceSize_;
170 
171  //- Number of times the element has been grown
172  label nGrowth_;
173 
174  //- Tolerance for the Ellipsoid of accuracy
175  static scalar tolerance_;
176 
177  //- Number of active species stored in the chemPoint
178  label nActiveSpecies_;
179 
180  //- Vectors that store the index conversion between the simplified
181  // and the complete chemical mechanism
182  List<label> simplifiedToCompleteIndex_;
183 
184  label timeTag_;
185  label lastTimeUsed_;
186 
187  bool toRemove_;
188 
189  label maxNumNewDim_;
190 
191  Switch printProportion_;
192 
193  //- Variable to store the number of retrieves the chemPoint
194  // will generate at each time step
195  label numRetrieve_;
196 
197  //- Variable to store the number of time steps the chempoint is allowed
198  // to still live according to the maxChPLifeTime_ parameter
199  label nLifeTime_;
200 
201  List<label> completeToSimplifiedIndex_;
202 
203  //- Number of equations in addition to the species eqs.
204  label nAdditionalEqns_;
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,
257  );
258 
259  //- Construct from another chemPoint and reference to a binary node
261  (
264  );
265 
266  //- Construct from another chemPoint
268  (
270  );
271 
272 
273  // Member functions
274 
275  //- Access to the TDACChemistryModel
277  {
278  return chemistry_;
279  }
281  inline label nGrowth()
282  {
283  return nGrowth_;
284  }
286  inline label& completeSpaceSize()
287  {
288  return completeSpaceSize_;
289  }
291  inline const scalarField& phi() const
292  {
293  return phi_;
294  }
296  inline const scalarField& Rphi() const
297  {
298  return Rphi_;
299  }
301  inline const scalarField& scaleFactor()
302  {
303  return scaleFactor_;
304  }
306  inline const scalar& tolerance()
307  {
308  return tolerance_;
309  }
311  static void changeTolerance(scalar newTol)
312  {
313  tolerance_ = newTol;
314  }
317  {
318  return node_;
319  }
321  inline const scalarSquareMatrix& A() const
322  {
323  return A_;
324  }
326  inline scalarSquareMatrix& A()
327  {
328  return A_;
329  }
331  inline const scalarSquareMatrix& LT() const
332  {
333  return LT_;
334  }
336  inline scalarSquareMatrix& LT()
337  {
338  return LT_;
339  }
341  inline label nActiveSpecies()
342  {
343  return nActiveSpecies_;
344  }
347  {
348  return completeToSimplifiedIndex_;
349  }
352  {
353  return simplifiedToCompleteIndex_;
354  }
355 
356  //- Increases the number of retrieves the chempoint has generated
357  void increaseNumRetrieve();
358 
359  //- Resets the number of retrieves at each time step
360  void resetNumRetrieve();
361 
362  //- Increases the "counter" of the chP life
363  void increaseNLifeTime();
364 
367  inline const label& timeTag()
368  {
369  return timeTag_;
370  }
372  inline label& lastTimeUsed()
373  {
374  return lastTimeUsed_;
375  }
377  inline bool& toRemove()
378  {
379  return toRemove_;
380  }
382  inline label& maxNumNewDim()
383  {
384  return maxNumNewDim_;
385  }
387  inline const label& numRetrieve()
388  {
389  return numRetrieve_;
390  }
392  inline const label& nLifeTime()
393  {
394  return nLifeTime_;
395  }
397  inline bool variableTimeStep() const
398  {
399  return chemistry_.variableTimeStep();
400  }
401 
402  // ISAT functions
403 
404  //- To RETRIEVE the mapping from the stored chemPoint phi, the query
405  // point phiq has to be in the EOA of phi.
406  // To test if phiq is in the ellipsoid:
407  // ||L^T.dphi|| <= 1
408  bool inEOA(const scalarField& phiq);
409 
410  //- More details about the minimum-volume ellipsoid covering an
411  // ellispoid E and a point p are found in [1].
412  // Here is the main steps to obtain the
413  // modified matrix L' describind the new ellipsoid.
414  // 1) calculate the point p' in the transformed space :
415  // p' = L^T.(p-phi)
416  // 2) compute the rank-one decomposition:
417  // G = I + gamma.p'.p'^T,
418  // with gamma = (1/|p'|-1)*1/|p'|^2
419  // 3) compute L':
420  // L'L'^T = (L.G)(L.G)^T,
421  // L'^T is then obtained by QR decomposition of (L.G)^T = G^T.L^T
422  // [1] Stephen B. Pope, "Algorithms for ellipsoids", FDA 08-01,
423  // Cornell University, 2008
424  bool grow(const scalarField& phiq);
425 
426  //- If phiq is not in the EOA, then the mapping is computed.
427  // But as the EOA is a conservative approximation of the region of
428  // accuracy surrounding the point phi, we could expand it by
429  // comparing thecomputed results with the one obtained by linear
430  // interpolation. The error eps is calculated:
431  // eps = ||B.(dR - dRl)||,
432  // with dR = Rphiq - Rphi, dRl = A.dphi and B the diagonal scale
433  // factor matrix.
434  // If eps <= tolerance, the EOA is too conservative and a GROW is
435  // performed,
436  // otherwise, the newly computed mapping is associated to the
437  // initial composition and added to the tree.
438  bool checkSolution
439  (
440  const scalarField& phiq,
441  const scalarField& Rphiq
442  );
443 };
444 
445 
446 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
447 
448 } // End namespace Foam
449 
450 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
451 
452 #ifdef NoRepository
453  #include "chemPointISAT.C"
454 #endif
455 
456 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
457 
458 #endif
459 
460 // ************************************************************************* //
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.
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()
binaryNode< CompType, ThermoType > *& node()
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
bool checkSolution(const scalarField &phiq, const scalarField &Rphiq)
If phiq is not in the EOA, then the mapping is computed.
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.
Definition: Switch.H:60
const label & timeTag()
const label & nLifeTime()
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
const scalarField & scaleFactor()
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
bool variableTimeStep() const
TDACChemistryModel< CompType, ThermoType > & chemistry()
Access to the TDACChemistryModel.
bool variableTimeStep() const
Return true if the time-step is variable and/or non-uniform.
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.