All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-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 \*---------------------------------------------------------------------------*/
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.mixture().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  const label li
64 )
65 {
66  scalarField c1(this->nSpecie_+2, 0.0);
67 
68  for(label i=0; i<this->nSpecie_; i++)
69  {
70  c1[i] = c[i];
71  }
72 
73  c1[this->nSpecie_] = T;
74  c1[this->nSpecie_+1] = p;
75 
76  // Compute the rAB matrix
77  RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
78  scalarField rABDen(this->nSpecie_,0.0);
79 
80  // Number of initialised rAB for each lines
81  Field<label> NbrABInit(this->nSpecie_,0);
82 
83  // Position of the initialised rAB, -1 when not initialised
84  RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
85 
86  // Index of the other species involved in the rABNum
87  RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
88 
89  forAll(this->chemistry_.reactions(), i)
90  {
91  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
92 
93  // For each reaction compute omegai
94  scalar omegaf, omegar;
95  const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
96 
97  // Then for each pair of species composing this reaction,
98  // compute the rAB matrix (separate the numerator and
99  // denominator)
100  DynamicList<scalar> wA(R.lhs().size()+R.rhs().size());
101  DynamicList<label> wAID(R.lhs().size()+R.rhs().size());
102 
103  forAll(R.lhs(), s)
104  {
105  label ss = R.lhs()[s].index;
106  scalar sl = -R.lhs()[s].stoichCoeff;
107  bool found(false);
108  forAll(wAID, id)
109  {
110  if (ss==wAID[id])
111  {
112  wA[id] += sl*omegai;
113  found = true;
114  break;
115  }
116  }
117  if (!found)
118  {
119  wA.append(sl*omegai);
120  wAID.append(ss);
121  }
122  }
123  forAll(R.rhs(), s)
124  {
125  label ss = R.rhs()[s].index;
126  scalar sl = R.rhs()[s].stoichCoeff;
127  bool found(false);
128  forAll(wAID, id)
129  {
130  if (ss==wAID[id])
131  {
132  wA[id] += sl*omegai;
133  found = true;
134  break;
135  }
136  }
137  if (!found)
138  {
139  wA.append(sl*omegai);
140  wAID.append(ss);
141  }
142  }
143 
144  // Now that all nuAi*wi are computed, without counting twice species
145  // present in both rhs and lhs, we can update the denominator and
146  // numerator for the rAB
147  wAID.shrink();
148  forAll(wAID, id)
149  {
150  label curID = wAID[id];
151 
152  // Absolute value of aggregated value
153  scalar curwA = ((wA[id]>=0) ? wA[id] : -wA[id]);
154 
155  List<bool> deltaBi(this->nSpecie_, false);
156  FIFOStack<label> usedIndex;
157  forAll(R.lhs(), j)
158  {
159  label sj = R.lhs()[j].index;
160  usedIndex.push(sj);
161  deltaBi[sj] = true;
162  }
163  forAll(R.rhs(), j)
164  {
165  label sj = R.rhs()[j].index;
166  usedIndex.push(sj);
167  deltaBi[sj] = true;
168  }
169 
170  // Disable for self reference (by definition rAA=0)
171  deltaBi[curID] = false;
172  while(!usedIndex.empty())
173  {
174  label curIndex = usedIndex.pop();
175 
176  if (deltaBi[curIndex])
177  {
178  // Disable to avoid counting it more than once
179  deltaBi[curIndex] = false;
180 
181  // Test if this rAB is not initialised
182  if (rABPos(curID, curIndex)==-1)
183  {
184  rABPos(curID, curIndex) = NbrABInit[curID];
185  NbrABInit[curID]++;
186  rABNum(curID, rABPos(curID, curIndex)) = curwA;
187  rABOtherSpec(curID, rABPos(curID, curIndex)) = curIndex;
188  }
189  else
190  {
191  rABNum(curID, rABPos(curID, curIndex)) += curwA;
192  }
193  }
194  }
195  if (rABDen[curID] == 0.0)
196  {
197  rABDen[curID] = curwA;
198  }
199  else
200  {
201  rABDen[curID] +=curwA;
202  }
203  }
204  }
205  // rii = 0.0 by definition
206 
207  label speciesNumber = 0;
208 
209  // Set all species to inactive and activate them according
210  // to rAB and initial set
211  for (label i=0; i<this->nSpecie_; i++)
212  {
213  this->activeSpecies_[i] = false;
214  }
215 
217 
218  // Initialise the list of active species with the search initiating set
219  // (SIS)
220  for (label i=0; i<searchInitSet_.size(); i++)
221  {
222  label q = searchInitSet_[i];
223  this->activeSpecies_[q] = true;
224  speciesNumber++;
225  Q.push(q);
226  }
227 
228  // Breadth first search with rAB
229  while (!Q.empty())
230  {
231  label u = Q.pop();
232  scalar Den = rABDen[u];
233 
234  if (Den > vSmall)
235  {
236  for (label v=0; v<NbrABInit[u]; v++)
237  {
238  label otherSpec = rABOtherSpec(u, v);
239  scalar rAB = rABNum(u, v)/Den;
240 
241  if (rAB > 1)
242  {
243  rAB = 1;
244  }
245 
246  // Include B only if rAB is above the tolerance and if the
247  // species was not searched before
248  if
249  (
250  rAB >= this->tolerance()
251  && !this->activeSpecies_[otherSpec]
252  )
253  {
254  Q.push(otherSpec);
255  this->activeSpecies_[otherSpec] = true;
256  speciesNumber++;
257  }
258  }
259  }
260  }
261 
262  // Put a flag on the reactions containing at least one removed species
263  forAll(this->chemistry_.reactions(), i)
264  {
265  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
266  this->chemistry_.reactionsDisabled()[i] = false;
267 
268  forAll(R.lhs(), s)
269  {
270  label ss = R.lhs()[s].index;
271 
272  // The species is inactive then the reaction is removed
273  if (!this->activeSpecies_[ss])
274  {
275  // Flag the reaction to disable it
276  this->chemistry_.reactionsDisabled()[i] = true;
277  break;
278  }
279  }
280 
281  // If the reaction has not been disabled yet
282  if (!this->chemistry_.reactionsDisabled()[i])
283  {
284  forAll(R.rhs(), s)
285  {
286  label ss = R.rhs()[s].index;
287  if (!this->activeSpecies_[ss])
288  {
289  this->chemistry_.reactionsDisabled()[i] = true;
290  break;
291  }
292  }
293  }
294  }
295 
296  this->NsSimp_ = speciesNumber;
297  this->chemistry_.simplifiedC().setSize(this->NsSimp_+2);
298  this->chemistry_.simplifiedToCompleteIndex().setSize(this->NsSimp_);
299 
300  label j = 0;
301  for (label i=0; i<this->nSpecie_; i++)
302  {
303  if (this->activeSpecies_[i])
304  {
305  this->chemistry_.simplifiedToCompleteIndex()[j] = i;
306  this->chemistry_.simplifiedC()[j] = c[i];
307  this->chemistry_.completeToSimplifiedIndex()[i] = j++;
308  if (!this->chemistry_.active(i))
309  {
310  this->chemistry_.setActive(i);
311  }
312  }
313  else
314  {
315  this->chemistry_.completeToSimplifiedIndex()[i] = -1;
316  }
317  }
318 
319  this->chemistry_.simplifiedC()[this->NsSimp_] = T;
320  this->chemistry_.simplifiedC()[this->NsSimp_+1] = p;
321  this->chemistry_.setNsDAC(this->NsSimp_);
322 
323  // Change temporary Ns in chemistryModel
324  // to make the function nEqns working
325  this->chemistry_.setNSpecie(this->NsSimp_);
326 }
327 
328 
329 // ************************************************************************* //
A FIFO stack based on a singly-linked list.
Definition: FIFOStack.H:51
Extends standardChemistryModel by adding the TDAC method.
A HashTable with keys but without contents.
Definition: HashSet.H:59
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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 List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: reactionI.H:36
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, const label li)
Reduce the mechanism.
Definition: DRG.C:59
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:55
void omega(const scalar p, const scalar T, const scalarField &c, const label li, scalarField &dcdt) const
Net reaction rate for individual species.
Definition: Reaction.C:290
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
basicChemistryModel & chemistry
const volScalarField & T
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
Definition: reactionI.H:42
DRG(const IOdictionary &dict, TDACChemistryModel< ThermoType > &chemistry)
Construct from components.
Definition: DRG.C:32
#define R(A, B, C, D, E, F, K, M)
void push(const T &a)
Push an element onto the stack.
Definition: FIFOStack.H:84
const multiComponentMixture< ThermoType > & mixture() const
Return reference to the mixture.
T pop()
Pop the bottom element off the stack.
Definition: FIFOStack.H:90
An abstract class for methods of chemical mechanism reduction.
volScalarField & p
bool found