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-2019 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 scalar p,
73  const scalar T,
74  const scalarField& c,
75  const label li
76 )
77 {
78  scalarField c1(this->nSpecie_+2, 0.0);
79 
80  for(label i=0; i<this->nSpecie_; i++)
81  {
82  c1[i] = c[i];
83  }
84 
85  c1[this->nSpecie_] = T;
86  c1[this->nSpecie_+1] = p;
87 
88  // Compute the rAB matrix
89  RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
90  scalarField rABDen(this->nSpecie_,0.0);
91 
92  // Number of initialized rAB for each lines
93  Field<label> NbrABInit(this->nSpecie_,0);
94 
95  // Position of the initialized rAB, -1 when not initialized
96  RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
97 
98  // Index of the other species involved in the rABNum
99  RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
100 
101  scalar pf, cf, pr, cr;
102  label lRef, rRef;
103  forAll(this->chemistry_.reactions(), i)
104  {
105  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
106 
107  // For each reaction compute omegai
108  scalar omegai = this->chemistry_.omega
109  (
110  R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef
111  );
112 
113 
114  // Then for each pair of species composing this reaction,
115  // compute the rAB matrix (separate the numerator and
116  // denominator)
117  DynamicList<scalar> wA(R.lhs().size()+R.rhs().size());
118  DynamicList<label> wAID(R.lhs().size()+R.rhs().size());
119 
120  forAll(R.lhs(), s)
121  {
122  label ss = R.lhs()[s].index;
123  scalar sl = -R.lhs()[s].stoichCoeff;
124  bool found(false);
125  forAll(wAID, id)
126  {
127  if (ss==wAID[id])
128  {
129  wA[id] += sl*omegai;
130  found = true;
131  break;
132  }
133  }
134  if (!found)
135  {
136  wA.append(sl*omegai);
137  wAID.append(ss);
138  }
139  }
140  forAll(R.rhs(), s)
141  {
142  label ss = R.rhs()[s].index;
143  scalar sl = R.rhs()[s].stoichCoeff;
144  bool found(false);
145  forAll(wAID, id)
146  {
147  if (ss==wAID[id])
148  {
149  wA[id] += sl*omegai;
150  found = true;
151  break;
152  }
153  }
154  if (!found)
155  {
156  wA.append(sl*omegai);
157  wAID.append(ss);
158  }
159  }
160 
161  // Now that all nuAi*wi are computed, without counting twice species
162  // present in both rhs and lhs, we can update the denominator and
163  // numerator for the rAB
164  wAID.shrink();
165  forAll(wAID, id)
166  {
167  label curID = wAID[id];
168 
169  // Absolute value of aggregated value
170  scalar curwA = ((wA[id]>=0) ? wA[id] : -wA[id]);
171 
172  List<bool> deltaBi(this->nSpecie_, false);
173  FIFOStack<label> usedIndex;
174  forAll(R.lhs(), j)
175  {
176  label sj = R.lhs()[j].index;
177  usedIndex.push(sj);
178  deltaBi[sj] = true;
179  }
180  forAll(R.rhs(), j)
181  {
182  label sj = R.rhs()[j].index;
183  usedIndex.push(sj);
184  deltaBi[sj] = true;
185  }
186 
187  // Disable for self reference (by definition rAA=0)
188  deltaBi[curID] = false;
189  while(!usedIndex.empty())
190  {
191  label curIndex = usedIndex.pop();
192 
193  if (deltaBi[curIndex])
194  {
195  // Disable to avoid counting it more than once
196  deltaBi[curIndex] = false;
197 
198  // Test if this rAB is not initialized
199  if (rABPos(curID, curIndex)==-1)
200  {
201  rABPos(curID, curIndex) = NbrABInit[curID];
202  NbrABInit[curID]++;
203  rABNum(curID, rABPos(curID, curIndex)) = curwA;
204  rABOtherSpec(curID, rABPos(curID, curIndex)) = curIndex;
205  }
206  else
207  {
208  rABNum(curID, rABPos(curID, curIndex)) += curwA;
209  }
210  }
211  }
212  if (rABDen[curID] == 0.0)
213  {
214  rABDen[curID] = curwA;
215  }
216  else
217  {
218  rABDen[curID] +=curwA;
219  }
220  }
221  }
222  // rii = 0.0 by definition
223 
224  label speciesNumber = 0;
225 
226  // Set all species to inactive and activate them according
227  // to rAB and initial set
228  for (label i=0; i<this->nSpecie_; i++)
229  {
230  this->activeSpecies_[i] = false;
231  }
232 
234 
235  // Initialize the list of active species with the search initiating set
236  // (SIS)
237  for (label i=0; i<searchInitSet_.size(); i++)
238  {
239  label q = searchInitSet_[i];
240  this->activeSpecies_[q] = true;
241  speciesNumber++;
242  Q.push(q);
243  }
244 
245  // Breadth first search with rAB
246  while (!Q.empty())
247  {
248  label u = Q.pop();
249  scalar Den = rABDen[u];
250 
251  if (Den > vSmall)
252  {
253  for (label v=0; v<NbrABInit[u]; v++)
254  {
255  label otherSpec = rABOtherSpec(u, v);
256  scalar rAB = rABNum(u, v)/Den;
257 
258  if (rAB > 1)
259  {
260  rAB = 1;
261  }
262 
263  // Include B only if rAB is above the tolerance and if the
264  // species was not searched before
265  if
266  (
267  rAB >= this->tolerance()
268  && !this->activeSpecies_[otherSpec]
269  )
270  {
271  Q.push(otherSpec);
272  this->activeSpecies_[otherSpec] = true;
273  speciesNumber++;
274  }
275  }
276  }
277  }
278 
279  // Put a flag on the reactions containing at least one removed species
280  forAll(this->chemistry_.reactions(), i)
281  {
282  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
283  this->chemistry_.reactionsDisabled()[i] = false;
284 
285  forAll(R.lhs(), s)
286  {
287  label ss = R.lhs()[s].index;
288 
289  // The species is inactive then the reaction is removed
290  if (!this->activeSpecies_[ss])
291  {
292  // Flag the reaction to disable it
293  this->chemistry_.reactionsDisabled()[i] = true;
294  break;
295  }
296  }
297 
298  // If the reaction has not been disabled yet
299  if (!this->chemistry_.reactionsDisabled()[i])
300  {
301  forAll(R.rhs(), s)
302  {
303  label ss = R.rhs()[s].index;
304  if (!this->activeSpecies_[ss])
305  {
306  this->chemistry_.reactionsDisabled()[i] = true;
307  break;
308  }
309  }
310  }
311  }
312 
313  this->NsSimp_ = speciesNumber;
314  this->chemistry_.simplifiedC().setSize(this->NsSimp_+2);
315  this->chemistry_.simplifiedToCompleteIndex().setSize(this->NsSimp_);
316 
317  label j = 0;
318  for (label i=0; i<this->nSpecie_; i++)
319  {
320  if (this->activeSpecies_[i])
321  {
322  this->chemistry_.simplifiedToCompleteIndex()[j] = i;
323  this->chemistry_.simplifiedC()[j] = c[i];
324  this->chemistry_.completeToSimplifiedIndex()[i] = j++;
325  if (!this->chemistry_.active(i))
326  {
327  this->chemistry_.setActive(i);
328  }
329  }
330  else
331  {
332  this->chemistry_.completeToSimplifiedIndex()[i] = -1;
333  }
334  }
335 
336  this->chemistry_.simplifiedC()[this->NsSimp_] = T;
337  this->chemistry_.simplifiedC()[this->NsSimp_+1] = p;
338  this->chemistry_.setNsDAC(this->NsSimp_);
339 
340  // Change temporary Ns in chemistryModel
341  // to make the function nEqns working
342  this->chemistry_.setNSpecie(this->NsSimp_);
343 }
344 
345 
346 // ************************************************************************* //
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:667
#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
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
const dimensionedScalar & c1
First radiation constant: default SI units: [W/m^2].
DRG(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
Definition: DRG.C:32
const PtrList< volScalarField > & Y()
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: ReactionI.H:57
BasicChemistryModel< rhoReactionThermo > & chemistry
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, const label li)
Reduce the mechanism.
Definition: DRG.C:71
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:934
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))
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
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