binaryNode.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 | Copyright (C) 2016-2017 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 "binaryNode.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class CompType, class ThermoType>
32 :
33  leafLeft_(nullptr),
34  leafRight_(nullptr),
35  nodeLeft_(nullptr),
36  nodeRight_(nullptr),
37  parent_(nullptr),
38  nAdditionalEqns_(0)
39 {}
40 
41 
42 template<class CompType, class ThermoType>
44 (
48 )
49 :
50  leafLeft_(elementLeft),
51  leafRight_(elementRight),
52  nodeLeft_(nullptr),
53  nodeRight_(nullptr),
54  parent_(parent),
55  v_(elementLeft->completeSpaceSize(), 0)
56 {
57  if (elementLeft->variableTimeStep())
58  {
59  nAdditionalEqns_ = 3;
60  }
61  else
62  {
63  nAdditionalEqns_ = 2;
64  }
65 
66  calcV(elementLeft, elementRight, v_);
67  a_ = calcA(elementLeft, elementRight);
68 }
69 
70 
71 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
72 
73 template<class CompType, class ThermoType>
75 (
78  scalarField& v
79 )
80 {
81  // LT is the transpose of the L matrix
82  scalarSquareMatrix& LT = elementLeft->LT();
83  bool mechReductionActive = elementLeft->chemistry().mechRed()->active();
84 
85  // Difference of composition in the full species domain
86  scalarField phiDif(elementRight->phi() - elementLeft->phi());
87  const scalarField& scaleFactor(elementLeft->scaleFactor());
88  scalar epsTol = elementLeft->tolerance();
89 
90  // v = LT.T()*LT*phiDif
91  for (label i=0; i<elementLeft->completeSpaceSize(); i++)
92  {
93  label si = i;
94  bool outOfIndexI = true;
95  if (mechReductionActive)
96  {
97  if (i<elementLeft->completeSpaceSize() - nAdditionalEqns_)
98  {
99  si = elementLeft->completeToSimplifiedIndex()[i];
100  outOfIndexI = (si == -1);
101  }
102  else // temperature and pressure
103  {
104  outOfIndexI = false;
105  const label dif =
106  i - (elementLeft->completeSpaceSize() - nAdditionalEqns_);
107  si = elementLeft->nActiveSpecies() + dif;
108  }
109  }
110  if (!mechReductionActive || (mechReductionActive && !(outOfIndexI)))
111  {
112  v[i] = 0;
113  for (label j=0; j<elementLeft->completeSpaceSize(); j++)
114  {
115  label sj = j;
116  bool outOfIndexJ = true;
117  if (mechReductionActive)
118  {
119  if (j < elementLeft->completeSpaceSize() - nAdditionalEqns_)
120  {
121  sj = elementLeft->completeToSimplifiedIndex()[j];
122  outOfIndexJ = (sj==-1);
123  }
124  else
125  {
126  outOfIndexJ = false;
127  const label dif =
128  j
129  - (
130  elementLeft->completeSpaceSize()
132  );
133  sj = elementLeft->nActiveSpecies() + dif;
134  }
135  }
136  if
137  (
138  !mechReductionActive
139  ||(mechReductionActive && !(outOfIndexJ))
140  )
141  {
142  // Since L is a lower triangular matrix k=0->min(i, j)
143  for (label k=0; k<=min(si, sj); k++)
144  {
145  v[i] += LT(k, si)*LT(k, sj)*phiDif[j];
146  }
147  }
148  }
149  }
150  else
151  {
152  // When it is an inactive species the diagonal element of LT is
153  // 1/(scaleFactor*epsTol)
154  v[i] = phiDif[i]/sqr(scaleFactor[i]*epsTol);
155  }
156  }
157 }
158 
159 
160 template<class CompType, class ThermoType>
162 (
165 )
166 {
167  scalarField phih((elementLeft->phi() + elementRight->phi())/2);
168  scalar a = 0;
169  forAll(phih, i)
170  {
171  a += v_[i]*phih[i];
172  }
173 
174  return a;
175 }
176 
177 
178 // ************************************************************************* //
scalar calcA(chemPointISAT< CompType, ThermoType > *elementLeft, chemPointISAT< CompType, ThermoType > *elementRight)
Compute a the product v^T.phih, with phih = (phi0 + phiq)/2.
Definition: binaryNode.C:162
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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 & LT() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Leaf of the binary tree. The chemPoint stores the composition &#39;phi&#39;, the mapping of this composition ...
List< label > & completeToSimplifiedIndex()
label k
Boltzmann constant.
scalarField v_
Definition: binaryNode.H:72
const scalar & a() const
Definition: binaryNode.H:179
binaryNode()
Construct null.
Definition: binaryNode.C:31
chemPointISAT< CompType, ThermoType > * leafLeft_
Element on the left.
Definition: binaryNode.H:55
binaryNode< CompType, ThermoType > * nodeLeft_
Node which follows on the left.
Definition: binaryNode.H:61
const scalarField & scaleFactor()
binaryNode< CompType, ThermoType > * nodeRight_
Node which follows on the right.
Definition: binaryNode.H:64
chemPointISAT< CompType, ThermoType > * leafRight_
Element on the right.
Definition: binaryNode.H:58
Node of the binary tree.
Definition: binaryNode.H:49
label nAdditionalEqns_
Number of equations in addition to the species eqs.
Definition: binaryNode.H:70
const scalarField & phi() const
bool variableTimeStep() const
TDACChemistryModel< CompType, ThermoType > & chemistry()
Access to the TDACChemistryModel.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void calcV(chemPointISAT< CompType, ThermoType > *&elementLeft, chemPointISAT< CompType, ThermoType > *&elementRight, scalarField &v)
Compute vector v:
Definition: binaryNode.C:75
label & completeSpaceSize()
binaryNode< CompType, ThermoType > * parent_
Parent node.
Definition: binaryNode.H:67