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