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-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 \*---------------------------------------------------------------------------*/
25 
26 #include "DRG.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class CompType, class ThermoType>
32 (
33  const IOdictionary& dict,
35 )
36 :
38  searchInitSet_(this->coeffsDict_.subDict("initialSet").size())
39 {
40  label j=0;
41  dictionary initSet = this->coeffsDict_.subDict("initialSet");
42  for (label i=0; i<chemistry.nSpecie(); i++)
43  {
44  if (initSet.found(chemistry.Y()[i].member()))
45  {
46  searchInitSet_[j++] = i;
47  }
48  }
49  if (j<searchInitSet_.size())
50  {
52  << searchInitSet_.size()-j
53  << " species in the initial set is not in the mechanism "
54  << initSet
55  << exit(FatalError);
56  }
57 }
58 
59 
60 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
61 
62 template<class CompType, class ThermoType>
64 {}
65 
66 
67 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
68 
69 template<class CompType, class ThermoType>
71 (
72  const scalarField &c,
73  const scalar T,
74  const scalar p
75 )
76 {
77  scalarField c1(this->nSpecie_+2, 0.0);
78 
79  for(label i=0; i<this->nSpecie_; i++)
80  {
81  c1[i] = c[i];
82  }
83 
84  c1[this->nSpecie_] = T;
85  c1[this->nSpecie_+1] = p;
86 
87  // Compute the rAB matrix
88  RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
89  scalarField rABDen(this->nSpecie_,0.0);
90 
91  // Number of initialized rAB for each lines
92  Field<label> NbrABInit(this->nSpecie_,0);
93 
94  // Position of the initialized rAB, -1 when not initialized
95  RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
96 
97  // Index of the other species involved in the rABNum
98  RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
99 
100  scalar pf, cf, pr, cr;
101  label lRef, rRef;
102  forAll(this->chemistry_.reactions(), i)
103  {
104  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
105  // For each reaction compute omegai
106  scalar omegai = this->chemistry_.omega
107  (
108  R, c1, T, p, pf, cf, lRef, pr, cr, rRef
109  );
110 
111 
112  // Then for each pair of species composing this reaction,
113  // compute the rAB matrix (separate the numerator and
114  // denominator)
115  DynamicList<scalar> wA(R.lhs().size()+R.rhs().size());
116  DynamicList<label> wAID(R.lhs().size()+R.rhs().size());
117 
118  forAll(R.lhs(), s)
119  {
120  label ss = R.lhs()[s].index;
121  scalar sl = -R.lhs()[s].stoichCoeff;
122  bool found(false);
123  forAll(wAID, id)
124  {
125  if (ss==wAID[id])
126  {
127  wA[id] += sl*omegai;
128  found = true;
129  break;
130  }
131  }
132  if (!found)
133  {
134  wA.append(sl*omegai);
135  wAID.append(ss);
136  }
137  }
138  forAll(R.rhs(), s)
139  {
140  label ss = R.rhs()[s].index;
141  scalar sl = R.rhs()[s].stoichCoeff;
142  bool found(false);
143  forAll(wAID, id)
144  {
145  if (ss==wAID[id])
146  {
147  wA[id] += sl*omegai;
148  found = true;
149  break;
150  }
151  }
152  if (!found)
153  {
154  wA.append(sl*omegai);
155  wAID.append(ss);
156  }
157  }
158 
159  // Now that all nuAi*wi are computed, without counting twice species
160  // present in both rhs and lhs, we can update the denominator and
161  // numerator for the rAB
162  wAID.shrink();
163  forAll(wAID, id)
164  {
165  label curID = wAID[id];
166 
167  // Absolute value of aggregated value
168  scalar curwA = ((wA[id]>=0) ? wA[id] : -wA[id]);
169 
170  List<bool> deltaBi(this->nSpecie_, false);
171  FIFOStack<label> usedIndex;
172  forAll(R.lhs(), j)
173  {
174  label sj = R.lhs()[j].index;
175  usedIndex.push(sj);
176  deltaBi[sj] = true;
177  }
178  forAll(R.rhs(), j)
179  {
180  label sj = R.rhs()[j].index;
181  usedIndex.push(sj);
182  deltaBi[sj] = true;
183  }
184 
185  // Disable for self reference (by definition rAA=0)
186  deltaBi[curID] = false;
187  while(!usedIndex.empty())
188  {
189  label curIndex = usedIndex.pop();
190 
191  if (deltaBi[curIndex])
192  {
193  // Disable to avoid counting it more than once
194  deltaBi[curIndex] = false;
195 
196  // Test if this rAB is not initialized
197  if (rABPos(curID, curIndex)==-1)
198  {
199  rABPos(curID, curIndex) = NbrABInit[curID];
200  NbrABInit[curID]++;
201  rABNum(curID, rABPos(curID, curIndex)) = curwA;
202  rABOtherSpec(curID, rABPos(curID, curIndex)) = curIndex;
203  }
204  else
205  {
206  rABNum(curID, rABPos(curID, curIndex)) += curwA;
207  }
208  }
209  }
210  if (rABDen[curID] == 0.0)
211  {
212  rABDen[curID] = curwA;
213  }
214  else
215  {
216  rABDen[curID] +=curwA;
217  }
218  }
219  }
220  // rii = 0.0 by definition
221 
222  label speciesNumber = 0;
223 
224  // Set all species to inactive and activate them according
225  // to rAB and initial set
226  for (label i=0; i<this->nSpecie_; i++)
227  {
228  this->activeSpecies_[i] = false;
229  }
230 
232 
233  // Initialize the list of active species with the search initiating set
234  // (SIS)
235  for (label i=0; i<searchInitSet_.size(); i++)
236  {
237  label q = searchInitSet_[i];
238  this->activeSpecies_[q] = true;
239  speciesNumber++;
240  Q.push(q);
241  }
242 
243  // Breadth first search with rAB
244  while (!Q.empty())
245  {
246  label u = Q.pop();
247  scalar Den = rABDen[u];
248 
249  if (Den > vSmall)
250  {
251  for (label v=0; v<NbrABInit[u]; v++)
252  {
253  label otherSpec = rABOtherSpec(u, v);
254  scalar rAB = rABNum(u, v)/Den;
255 
256  if (rAB > 1)
257  {
258  rAB = 1;
259  }
260 
261  // Include B only if rAB is above the tolerance and if the
262  // species was not searched before
263  if
264  (
265  rAB >= this->tolerance()
266  && !this->activeSpecies_[otherSpec]
267  )
268  {
269  Q.push(otherSpec);
270  this->activeSpecies_[otherSpec] = true;
271  speciesNumber++;
272  }
273  }
274  }
275  }
276 
277  // Put a flag on the reactions containing at least one removed species
278  forAll(this->chemistry_.reactions(), i)
279  {
280  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
281  this->chemistry_.reactionsDisabled()[i] = false;
282 
283  forAll(R.lhs(), s)
284  {
285  label ss = R.lhs()[s].index;
286 
287  // The species is inactive then the reaction is removed
288  if (!this->activeSpecies_[ss])
289  {
290  // Flag the reaction to disable it
291  this->chemistry_.reactionsDisabled()[i] = true;
292  break;
293  }
294  }
295 
296  // If the reaction has not been disabled yet
297  if (!this->chemistry_.reactionsDisabled()[i])
298  {
299  forAll(R.rhs(), s)
300  {
301  label ss = R.rhs()[s].index;
302  if (!this->activeSpecies_[ss])
303  {
304  this->chemistry_.reactionsDisabled()[i] = true;
305  break;
306  }
307  }
308  }
309  }
310 
311  this->NsSimp_ = speciesNumber;
312  this->chemistry_.simplifiedC().setSize(this->NsSimp_+2);
313  this->chemistry_.simplifiedToCompleteIndex().setSize(this->NsSimp_);
314 
315  label j = 0;
316  for (label i=0; i<this->nSpecie_; i++)
317  {
318  if (this->activeSpecies_[i])
319  {
320  this->chemistry_.simplifiedToCompleteIndex()[j] = i;
321  this->chemistry_.simplifiedC()[j] = c[i];
322  this->chemistry_.completeToSimplifiedIndex()[i] = j++;
323  if (!this->chemistry_.active(i))
324  {
325  this->chemistry_.setActive(i);
326  }
327  }
328  else
329  {
330  this->chemistry_.completeToSimplifiedIndex()[i] = -1;
331  }
332  }
333 
334  this->chemistry_.simplifiedC()[this->NsSimp_] = T;
335  this->chemistry_.simplifiedC()[this->NsSimp_+1] = p;
336  this->chemistry_.setNsDAC(this->NsSimp_);
337 
338  // Change temporary Ns in chemistryModel
339  // to make the function nEqns working
340  this->chemistry_.setNSpecie(this->NsSimp_);
341 }
342 
343 
344 // ************************************************************************* //
A FIFO stack based on a singly-linked list.
Definition: FIFOStack.H:51
virtual label nSpecie() const
The number of species.
dictionary dict
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:438
#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
PtrList< volScalarField > & Y()
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
DRG(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
Definition: DRG.C:32
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: ReactionI.H:57
BasicChemistryModel< rhoReactionThermo > & chemistry
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:56
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:52
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:699
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))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
Definition: DRG.C:71
const volScalarField & T
#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 dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
T pop()
Pop the bottom element off the stack.
Definition: FIFOStack.H:90
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
Definition: ReactionI.H:64
An abstract class for methods of chemical mechanism reduction.
volScalarField & p
bool found