DRG.C
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-2023 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 "DRG.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class ThermoType>
32 (
33  const IOdictionary& dict,
35 )
36 :
38  searchInitSet_()
39 {
40  const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
41  forAllConstIter(wordHashSet, initSet, iter)
42  {
43  searchInitSet_.append(chemistry.thermo().species()[iter.key()]);
44  }
45 }
46 
47 
48 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
49 
50 template<class ThermoType>
52 {}
53 
54 
55 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
56 
57 template<class ThermoType>
59 (
60  const scalar p,
61  const scalar T,
62  const scalarField& c,
63  List<label>& ctos,
64  DynamicList<label>& stoc,
65  const label li
66 )
67 {
69 
70  scalarField c1(this->nSpecie()+2, 0.0);
71 
72  for(label i=0; i<this->nSpecie(); i++)
73  {
74  c1[i] = c[i];
75  }
76 
77  c1[this->nSpecie()] = T;
78  c1[this->nSpecie()+1] = p;
79 
80  // Compute the rAB matrix
81  RectangularMatrix<scalar> rABNum(this->nSpecie(),this->nSpecie(),0.0);
82  scalarField rABDen(this->nSpecie(),0.0);
83 
84  // Number of initialised rAB for each lines
85  Field<label> NbrABInit(this->nSpecie(),0);
86 
87  // Position of the initialised rAB, -1 when not initialised
88  RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
89 
90  // Index of the other species involved in the rABNum
91  RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
92 
93  forAll(this->chemistry_.reactions(), i)
94  {
95  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
96 
97  // For each reaction compute omegai
98  scalar omegaf, omegar;
99  const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
100 
101  // Then for each pair of species composing this reaction,
102  // compute the rAB matrix (separate the numerator and
103  // denominator)
104  DynamicList<scalar> wA(R.lhs().size()+R.rhs().size());
105  DynamicList<label> wAID(R.lhs().size()+R.rhs().size());
106 
107  forAll(R.lhs(), s)
108  {
109  label ss = R.lhs()[s].index;
110  scalar sl = -R.lhs()[s].stoichCoeff;
111  bool found(false);
112  forAll(wAID, id)
113  {
114  if (ss==wAID[id])
115  {
116  wA[id] += sl*omegai;
117  found = true;
118  break;
119  }
120  }
121  if (!found)
122  {
123  wA.append(sl*omegai);
124  wAID.append(ss);
125  }
126  }
127  forAll(R.rhs(), s)
128  {
129  label ss = R.rhs()[s].index;
130  scalar sl = R.rhs()[s].stoichCoeff;
131  bool found(false);
132  forAll(wAID, id)
133  {
134  if (ss==wAID[id])
135  {
136  wA[id] += sl*omegai;
137  found = true;
138  break;
139  }
140  }
141  if (!found)
142  {
143  wA.append(sl*omegai);
144  wAID.append(ss);
145  }
146  }
147 
148  // Now that all nuAi*wi are computed, without counting twice species
149  // present in both rhs and lhs, we can update the denominator and
150  // numerator for the rAB
151  wAID.shrink();
152  forAll(wAID, id)
153  {
154  label curID = wAID[id];
155 
156  // Absolute value of aggregated value
157  scalar curwA = ((wA[id]>=0) ? wA[id] : -wA[id]);
158 
159  List<bool> deltaBi(this->nSpecie(), false);
160  FIFOStack<label> usedIndex;
161  forAll(R.lhs(), j)
162  {
163  label sj = R.lhs()[j].index;
164  usedIndex.push(sj);
165  deltaBi[sj] = true;
166  }
167  forAll(R.rhs(), j)
168  {
169  label sj = R.rhs()[j].index;
170  usedIndex.push(sj);
171  deltaBi[sj] = true;
172  }
173 
174  // Disable for self reference (by definition rAA=0)
175  deltaBi[curID] = false;
176  while(!usedIndex.empty())
177  {
178  label curIndex = usedIndex.pop();
179 
180  if (deltaBi[curIndex])
181  {
182  // Disable to avoid counting it more than once
183  deltaBi[curIndex] = false;
184 
185  // Test if this rAB is not initialised
186  if (rABPos(curID, curIndex)==-1)
187  {
188  rABPos(curID, curIndex) = NbrABInit[curID];
189  NbrABInit[curID]++;
190  rABNum(curID, rABPos(curID, curIndex)) = curwA;
191  rABOtherSpec(curID, rABPos(curID, curIndex)) = curIndex;
192  }
193  else
194  {
195  rABNum(curID, rABPos(curID, curIndex)) += curwA;
196  }
197  }
198  }
199  if (rABDen[curID] == 0.0)
200  {
201  rABDen[curID] = curwA;
202  }
203  else
204  {
205  rABDen[curID] +=curwA;
206  }
207  }
208  }
209  // rii = 0.0 by definition
210 
211  // Set all species to inactive and activate them according
212  // to rAB and initial set
213  for (label i=0; i<this->nSpecie(); i++)
214  {
215  this->activeSpecies_[i] = false;
216  }
217 
219 
220  // Initialise the list of active species with the search initiating set
221  // (SIS)
222  for (label i=0; i<searchInitSet_.size(); i++)
223  {
224  label q = searchInitSet_[i];
225  this->activeSpecies_[q] = true;
226  Q.push(q);
227  }
228 
229  // Breadth first search with rAB
230  while (!Q.empty())
231  {
232  label u = Q.pop();
233  scalar Den = rABDen[u];
234 
235  if (Den > vSmall)
236  {
237  for (label v=0; v<NbrABInit[u]; v++)
238  {
239  label otherSpec = rABOtherSpec(u, v);
240  scalar rAB = rABNum(u, v)/Den;
241 
242  if (rAB > 1)
243  {
244  rAB = 1;
245  }
246 
247  // Include B only if rAB is above the tolerance and if the
248  // species was not searched before
249  if
250  (
251  rAB >= this->tolerance()
252  && !this->activeSpecies_[otherSpec]
253  )
254  {
255  Q.push(otherSpec);
256  this->activeSpecies_[otherSpec] = true;
257  }
258  }
259  }
260  }
261 
263 }
264 
265 
266 // ************************************************************************* //
bool found
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:252
A FIFO stack based on a singly-linked list.
Definition: FIFOStack.H:54
void push(const T &a)
Push an element onto the stack.
Definition: FIFOStack.H:84
T pop()
Pop the bottom element off the stack.
Definition: FIFOStack.H:90
A HashTable with keys but without contents.
Definition: HashSet.H:62
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
Simple extension of ThermoType to handle reaction kinetics in addition to the equilibrium thermodynam...
Definition: Reaction.H:72
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
An abstract class for methods of chemical mechanism reduction.
void initReduceMechanism()
Protected Member Functions.
const dictionary coeffsDict_
Dictionary that store the algorithm data.
void endReduceMechanism(List< label > &ctos, DynamicList< label > &stoc)
End reduction of the mechanism.
DRG(const IOdictionary &dict, chemistryModel< ThermoType > &chemistry)
Construct from components.
Definition: DRG.C:32
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, List< label > &ctos, DynamicList< label > &stoc, const label li)
Reduce the mechanism.
Definition: DRG.C:59
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:710
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
const dimensionedScalar c
Speed of light in a vacuum.
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
static scalar R(const scalar a, const scalar x)
Definition: invIncGamma.C:102
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
const label nSpecie
dictionary dict
volScalarField & p
basicChemistryModel & chemistry