All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PFA.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 "PFA.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 
42  dictionary initSet = this->coeffsDict_.subDict("initialSet");
43 
44  for (label i=0; i<chemistry.nSpecie(); i++)
45  {
46  if (initSet.found(chemistry.Y()[i].member()))
47  {
48  searchInitSet_[j++] = i;
49  }
50  }
51 
52  if (j<searchInitSet_.size())
53  {
55  << searchInitSet_.size()-j
56  << " species in the initial set is not in the mechanism "
57  << initSet
58  << abort(FatalError);
59  }
60 }
61 
62 
63 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
64 
65 template<class CompType, class ThermoType>
67 {}
68 
69 
70 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
71 
72 template<class CompType, class ThermoType>
74 (
75  const scalarField &c,
76  const scalar T,
77  const scalar p
78 )
79 {
80  scalarField& completeC(this->chemistry_.completeC());
81  scalarField c1(this->chemistry_.nEqns(), 0.0);
82 
83  for (label i=0; i<this->nSpecie_; i++)
84  {
85  c1[i] = c[i];
86  completeC[i] = c[i];
87  }
88 
89  c1[this->nSpecie_] = T;
90  c1[this->nSpecie_+1] = p;
91 
92  // Compute the rAB matrix
93  RectangularMatrix<scalar> PAB(this->nSpecie_,this->nSpecie_,0.0);
94  RectangularMatrix<scalar> CAB(this->nSpecie_,this->nSpecie_,0.0);
95  scalarField PA(this->nSpecie_,0.0);
96  scalarField CA(this->nSpecie_,0.0);
97 
98  // Number of initialized rAB for each lines
99  Field<label> NbrABInit(this->nSpecie_,0);
100  // Position of the initialized rAB, -1 when not initialized
101  RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
102  // Index of the other species involved in the rABNum
103  RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
104 
105  scalar pf, cf, pr, cr;
106  label lRef, rRef;
107  forAll(this->chemistry_.reactions(), i)
108  {
109  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
110  // for each reaction compute omegai
111  scalar omegai = this->chemistry_.omega
112  (
113  R, c1, T, p, pf, cf, lRef, pr, cr, rRef
114  );
115 
116  // then for each pair of species composing this reaction,
117  // compute the rAB matrix (separate the numerator and
118  // denominator)
119 
120  DynamicList<scalar> wA(R.lhs().size()+R.rhs().size());
121  DynamicList<label> wAID(R.lhs().size()+R.rhs().size());
122 
123  forAll(R.lhs(), s)// compute rAB for all species in the left hand side
124  {
125  label ss = R.lhs()[s].index;
126  scalar sl = -R.lhs()[s].stoichCoeff; // vAi = v''-v' => here -v'
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  forAll(R.rhs(), s) // compute rAB for all species in the right hand side
144  {
145  label ss = R.rhs()[s].index;
146  scalar sl = R.rhs()[s].stoichCoeff; // vAi = v''-v' => here v''
147  bool found(false);
148  forAll(wAID, id)
149  {
150  if (ss==wAID[id])
151  {
152  wA[id] += sl*omegai;
153  found = true;
154  break;
155  }
156  }
157  if (!found)
158  {
159  wA.append(sl*omegai);
160  wAID.append(ss);
161  }
162  }
163 
164  wAID.shrink();
165  forAll(wAID, id)
166  {
167  label curID = wAID[id];
168  scalar curwA = wA[id];
169  List<bool> deltaBi(this->nSpecie_, false);
170  FIFOStack<label> usedIndex;
171  forAll(R.lhs(),j)
172  {
173  label sj = R.lhs()[j].index;
174  usedIndex.push(sj);
175  deltaBi[sj] = true;
176  }
177  forAll(R.rhs(),j)
178  {
179  label sj = R.rhs()[j].index;
180  usedIndex.push(sj);
181  deltaBi[sj] = true;
182  }
183 
184  deltaBi[curID] = false;
185 
186  while(!usedIndex.empty())
187  {
188  label curIndex = usedIndex.pop();
189 
190  if (deltaBi[curIndex])
191  {
192  deltaBi[curIndex] = false;
193  if (rABPos(curID, curIndex)==-1)
194  {
195  rABPos(curID, curIndex) = NbrABInit[curID];
196  rABOtherSpec(curID, NbrABInit[curID]) = curIndex;
197  if (curwA > 0.0)
198  {
199  PAB(curID, NbrABInit[curID]) = curwA;
200  }
201  else
202  {
203  CAB(curID, NbrABInit[curID]) = -curwA;
204  }
205  NbrABInit[curID]++;
206  }
207  else
208  {
209  if (curwA > 0.0)
210  {
211  PAB(curID, rABPos(curID, curIndex)) += curwA;
212  }
213  else
214  {
215  CAB(curID, rABPos(curID, curIndex)) += -curwA;
216  }
217  }
218  }
219  }
220  // Now that every species of the reactions has been visited, we can
221  // compute the production and consumption rate. It avoids getting
222  // wrong results when species are present in both lhs and rhs
223  if (curwA > 0.0)
224  {
225  if (PA[curID] == 0.0)
226  {
227  PA[curID] = curwA;
228  }
229  else
230  {
231  PA[curID] += curwA;
232  }
233  }
234  else
235  {
236  if (CA[curID] == 0.0)
237  {
238  CA[curID] = -curwA;
239  }
240  else
241  {
242  CA[curID] += -curwA;
243  }
244  }
245  }
246  }
247  // rii = 0.0 by definition
248 
249  // Compute second generation link strength
250  // For all species A look at all rAri of the connected species ri and
251  // compute rriB with all the connected species of ri, B different of A. If
252  // a new species is connected, add it to the list of connected species. It
253  // is a connection of second generation and it will be aggregated in the
254  // final step to evaluate the total connection strength (or path flux).
255  // Compute rsecond=rAri*rriB with A!=ri!=B
256  RectangularMatrix<scalar> PAB2nd(this->nSpecie_,this->nSpecie_,0.0);
257  RectangularMatrix<scalar> CAB2nd(this->nSpecie_,this->nSpecie_,0.0);
258 
259  // Number of initialized rAB for each lines
260  Field<label> NbrABInit2nd(this->nSpecie_, 0);
261 
262  // Position of the initialized rAB, -1 when not initialized
263  RectangularMatrix<label> rABPos2nd(this->nSpecie_, this->nSpecie_, -1);
264 
265  // Index of the other species involved in the rABNum
266  RectangularMatrix<label> rABOtherSpec2nd
267  (
268  this->nSpecie_, this->nSpecie_, -1
269  );
270 
271  forAll(NbrABInit, A)
272  {
273  for (int i=0; i<NbrABInit[A]; i++)
274  {
275  label ri = rABOtherSpec(A, i);
276  scalar maxPACA = max(PA[ri],CA[ri]);
277  if (maxPACA > vSmall)
278  {
279  for (int j=0; j<NbrABInit[ri]; j++)
280  {
281  label B = rABOtherSpec(ri, j);
282  if (B != A) // if B!=A and also !=ri by definition
283  {
284  if (rABPos2nd(A, B)==-1)
285  {
286  rABPos2nd(A, B) = NbrABInit2nd[A]++;
287  rABOtherSpec2nd(A, rABPos2nd(A, B)) = B;
288  PAB2nd(A, rABPos2nd(A, B)) =
289  PAB(A, i)*PAB(ri, j)/maxPACA;
290  CAB2nd(A, rABPos2nd(A, B)) =
291  CAB(A, i)*CAB(ri, j)/maxPACA;
292  }
293  else
294  {
295  PAB2nd(A, rABPos2nd(A, B)) +=
296  PAB(A, i)*PAB(ri, j)/maxPACA;
297  CAB2nd(A, rABPos2nd(A, B)) +=
298  CAB(A, i)*CAB(ri, j)/maxPACA;
299  }
300  }
301  }
302  }
303  }
304  }
305 
306  // Using the rAB matrix (numerator and denominator separated)
307  label speciesNumber = 0;
308 
309  // set all species to inactive and activate them according
310  // to rAB and initial set
311  for (label i=0; i<this->nSpecie_; i++)
312  {
313  this->activeSpecies_[i] = false;
314  }
315 
316  // Initialize the FIFOStack for search set
317  const labelList& SIS(this->searchInitSet_);
319 
320  for (label i=0; i<SIS.size(); i++)
321  {
322  label q = SIS[i];
323  this->activeSpecies_[q] = true;
324  speciesNumber++;
325  Q.push(q);
326  }
327 
328  // Execute the main loop for R-value
329  while (!Q.empty())
330  {
331  label u = Q.pop();
332  scalar Den = max(PA[u],CA[u]);
333 
334  if (Den!=0.0)
335  {
336  // first generation
337  for (label v=0; v<NbrABInit[u]; v++)
338  {
339  label otherSpec = rABOtherSpec(u, v);
340  scalar rAB = (PAB(u, v)+CAB(u, v))/Den; // first generation
341  label id2nd = rABPos2nd(u, otherSpec);
342  if (id2nd !=-1)// if there is a second generation link
343  {
344  rAB += (PAB2nd(u, id2nd)+CAB2nd(u, id2nd))/Den;
345  }
346  // the link is stronger than the user-defined tolerance
347  if
348  (
349  rAB >= this->tolerance()
350  && !this->activeSpecies_[otherSpec]
351  )
352  {
353  Q.push(otherSpec);
354  this->activeSpecies_[otherSpec] = true;
355  speciesNumber++;
356  }
357 
358  }
359  // second generation link only (for those without first link)
360  for (label v=0; v<NbrABInit2nd[u]; v++)
361  {
362  label otherSpec = rABOtherSpec2nd(u, v);
363  scalar rAB = (PAB2nd(u, v)+CAB2nd(u, v))/Den;
364  // the link is stronger than the user-defined tolerance
365  if
366  (
367  rAB >= this->tolerance()
368  && !this->activeSpecies_[otherSpec]
369  )
370  {
371  Q.push(otherSpec);
372  this->activeSpecies_[otherSpec] = true;
373  speciesNumber++;
374  }
375  }
376  }
377  }
378 
379  // Put a flag on the reactions containing at least one removed species
380  forAll(this->chemistry_.reactions(), i)
381  {
382  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
383  this->chemistry_.reactionsDisabled()[i] = false;
384 
385  forAll(R.lhs(), s)
386  {
387  label ss = R.lhs()[s].index;
388  if (!this->activeSpecies_[ss])
389  {
390  this->chemistry_.reactionsDisabled()[i] = true;
391  break;
392  }
393  }
394  if (!this->chemistry_.reactionsDisabled()[i])
395  {
396  forAll(R.rhs(), s)
397  {
398  label ss = R.rhs()[s].index;
399  if (!this->activeSpecies_[ss])
400  {
401  this->chemistry_.reactionsDisabled()[i] = true;
402  break;
403  }
404  }
405  }
406  }
407 
408  this->NsSimp_ = speciesNumber;
409  scalarField& simplifiedC(this->chemistry_.simplifiedC());
410  simplifiedC.setSize(this->NsSimp_+2);
411  DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
412  s2c.setSize(this->NsSimp_);
413  Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
414 
415  label j = 0;
416  for (label i=0; i<this->nSpecie_; i++)
417  {
418  if (this->activeSpecies_[i])
419  {
420  s2c[j] = i;
421  simplifiedC[j] = c[i];
422  c2s[i] = j++;
423  if (!this->chemistry_.active(i))
424  {
425  this->chemistry_.setActive(i);
426  }
427  }
428  else
429  {
430  c2s[i] = -1;
431  }
432  }
433  simplifiedC[this->NsSimp_] = T;
434  simplifiedC[this->NsSimp_+1] = p;
435  this->chemistry_.setNsDAC(this->NsSimp_);
436  // change temporary Ns in chemistryModel
437  // to make the function nEqns working
438  this->chemistry_.setNSpecie(this->NsSimp_);
439 }
440 
441 
442 // ************************************************************************* //
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()
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
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))
void setSize(const label)
Alter the addressed list size.
Definition: DynamicListI.H:175
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const volScalarField & T
virtual ~PFA()
Destructor.
Definition: PFA.C:66
void setSize(const label)
Reset size of List.
Definition: List.C:281
#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].
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
Definition: PFA.C:74
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
PFA(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
Definition: PFA.C:32
bool found