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-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 "PFA.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  List<label>& ctos,
64  DynamicList<label>& stoc,
65  const label li
66 )
67 {
69 
70  scalarField c1(this->chemistry_.nEqns(), 0.0);
71 
72  for (label i=0; i<this->nSpecie(); i++)
73  {
74  c1[i] = c[i];
75  }
76 
77  c1[this->nSpecie()] = T;
78  c1[this->nSpecie()+1] = p;
79 
80  // Compute the rAB matrix
81  RectangularMatrix<scalar> PAB(this->nSpecie(),this->nSpecie(),0.0);
82  RectangularMatrix<scalar> CAB(this->nSpecie(),this->nSpecie(),0.0);
83  scalarField PA(this->nSpecie(),0.0);
84  scalarField CA(this->nSpecie(),0.0);
85 
86  // Number of initialised rAB for each lines
87  Field<label> NbrABInit(this->nSpecie(),0);
88  // Position of the initialised rAB, -1 when not initialised
89  RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
90  // Index of the other species involved in the rABNum
91  RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
92 
93  forAll(this->chemistry_.reactions(), i)
94  {
95  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
96 
97  // for each reaction compute omegai
98  scalar omegaf, omegar;
99  const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
100 
101  // then for each pair of species composing this reaction,
102  // compute the rAB matrix (separate the numerator and
103  // denominator)
104 
105  DynamicList<scalar> wA(R.lhs().size()+R.rhs().size());
106  DynamicList<label> wAID(R.lhs().size()+R.rhs().size());
107 
108  forAll(R.lhs(), s)// compute rAB for all species in the left hand side
109  {
110  label ss = R.lhs()[s].index;
111  scalar sl = -R.lhs()[s].stoichCoeff; // vAi = v''-v' => here -v'
112  bool found(false);
113  forAll(wAID, id)
114  {
115  if (ss==wAID[id])
116  {
117  wA[id] += sl*omegai;
118  found = true;
119  break;
120  }
121  }
122  if (!found)
123  {
124  wA.append(sl*omegai);
125  wAID.append(ss);
126  }
127  }
128  forAll(R.rhs(), s) // compute rAB for all species in the right hand side
129  {
130  label ss = R.rhs()[s].index;
131  scalar sl = R.rhs()[s].stoichCoeff; // vAi = v''-v' => here v''
132  bool found(false);
133  forAll(wAID, id)
134  {
135  if (ss==wAID[id])
136  {
137  wA[id] += sl*omegai;
138  found = true;
139  break;
140  }
141  }
142  if (!found)
143  {
144  wA.append(sl*omegai);
145  wAID.append(ss);
146  }
147  }
148 
149  wAID.shrink();
150  forAll(wAID, id)
151  {
152  label curID = wAID[id];
153  scalar curwA = wA[id];
154  List<bool> deltaBi(this->nSpecie(), false);
155  FIFOStack<label> usedIndex;
156  forAll(R.lhs(),j)
157  {
158  label sj = R.lhs()[j].index;
159  usedIndex.push(sj);
160  deltaBi[sj] = true;
161  }
162  forAll(R.rhs(),j)
163  {
164  label sj = R.rhs()[j].index;
165  usedIndex.push(sj);
166  deltaBi[sj] = true;
167  }
168 
169  deltaBi[curID] = false;
170 
171  while(!usedIndex.empty())
172  {
173  label curIndex = usedIndex.pop();
174 
175  if (deltaBi[curIndex])
176  {
177  deltaBi[curIndex] = false;
178  if (rABPos(curID, curIndex)==-1)
179  {
180  rABPos(curID, curIndex) = NbrABInit[curID];
181  rABOtherSpec(curID, NbrABInit[curID]) = curIndex;
182  if (curwA > 0.0)
183  {
184  PAB(curID, NbrABInit[curID]) = curwA;
185  }
186  else
187  {
188  CAB(curID, NbrABInit[curID]) = -curwA;
189  }
190  NbrABInit[curID]++;
191  }
192  else
193  {
194  if (curwA > 0.0)
195  {
196  PAB(curID, rABPos(curID, curIndex)) += curwA;
197  }
198  else
199  {
200  CAB(curID, rABPos(curID, curIndex)) += -curwA;
201  }
202  }
203  }
204  }
205  // Now that every species of the reactions has been visited, we can
206  // compute the production and consumption rate. It avoids getting
207  // wrong results when species are present in both lhs and rhs
208  if (curwA > 0.0)
209  {
210  if (PA[curID] == 0.0)
211  {
212  PA[curID] = curwA;
213  }
214  else
215  {
216  PA[curID] += curwA;
217  }
218  }
219  else
220  {
221  if (CA[curID] == 0.0)
222  {
223  CA[curID] = -curwA;
224  }
225  else
226  {
227  CA[curID] += -curwA;
228  }
229  }
230  }
231  }
232  // rii = 0.0 by definition
233 
234  // Compute second generation link strength
235  // For all species A look at all rAri of the connected species ri and
236  // compute rriB with all the connected species of ri, B different of A. If
237  // a new species is connected, add it to the list of connected species. It
238  // is a connection of second generation and it will be aggregated in the
239  // final step to evaluate the total connection strength (or path flux).
240  // Compute rsecond=rAri*rriB with A!=ri!=B
241  RectangularMatrix<scalar> PAB2nd(this->nSpecie(),this->nSpecie(),0.0);
242  RectangularMatrix<scalar> CAB2nd(this->nSpecie(),this->nSpecie(),0.0);
243 
244  // Number of initialised rAB for each lines
245  Field<label> NbrABInit2nd(this->nSpecie(), 0);
246 
247  // Position of the initialised rAB, -1 when not initialised
248  RectangularMatrix<label> rABPos2nd(this->nSpecie(), this->nSpecie(), -1);
249 
250  // Index of the other species involved in the rABNum
251  RectangularMatrix<label> rABOtherSpec2nd
252  (
253  this->nSpecie(), this->nSpecie(), -1
254  );
255 
256  forAll(NbrABInit, A)
257  {
258  for (int i=0; i<NbrABInit[A]; i++)
259  {
260  label ri = rABOtherSpec(A, i);
261  scalar maxPACA = max(PA[ri],CA[ri]);
262  if (maxPACA > vSmall)
263  {
264  for (int j=0; j<NbrABInit[ri]; j++)
265  {
266  label B = rABOtherSpec(ri, j);
267  if (B != A) // if B!=A and also !=ri by definition
268  {
269  if (rABPos2nd(A, B)==-1)
270  {
271  rABPos2nd(A, B) = NbrABInit2nd[A]++;
272  rABOtherSpec2nd(A, rABPos2nd(A, B)) = B;
273  PAB2nd(A, rABPos2nd(A, B)) =
274  PAB(A, i)*PAB(ri, j)/maxPACA;
275  CAB2nd(A, rABPos2nd(A, B)) =
276  CAB(A, i)*CAB(ri, j)/maxPACA;
277  }
278  else
279  {
280  PAB2nd(A, rABPos2nd(A, B)) +=
281  PAB(A, i)*PAB(ri, j)/maxPACA;
282  CAB2nd(A, rABPos2nd(A, B)) +=
283  CAB(A, i)*CAB(ri, j)/maxPACA;
284  }
285  }
286  }
287  }
288  }
289  }
290 
291  // Using the rAB matrix (numerator and denominator separated)
292  // set all species to inactive and activate them according
293  // to rAB and initial set
294  for (label i=0; i<this->nSpecie(); i++)
295  {
296  this->activeSpecies_[i] = false;
297  }
298 
299  // Initialise the FIFOStack for search set
300  const labelList& SIS(this->searchInitSet_);
302 
303  for (label i=0; i<SIS.size(); i++)
304  {
305  label q = SIS[i];
306  this->activeSpecies_[q] = true;
307  Q.push(q);
308  }
309 
310  // Execute the main loop for R-value
311  while (!Q.empty())
312  {
313  label u = Q.pop();
314  scalar Den = max(PA[u],CA[u]);
315 
316  if (Den!=0.0)
317  {
318  // first generation
319  for (label v=0; v<NbrABInit[u]; v++)
320  {
321  label otherSpec = rABOtherSpec(u, v);
322  scalar rAB = (PAB(u, v)+CAB(u, v))/Den; // first generation
323  label id2nd = rABPos2nd(u, otherSpec);
324  if (id2nd !=-1)// if there is a second generation link
325  {
326  rAB += (PAB2nd(u, id2nd)+CAB2nd(u, id2nd))/Den;
327  }
328  // the link is stronger than the user-defined tolerance
329  if
330  (
331  rAB >= this->tolerance()
332  && !this->activeSpecies_[otherSpec]
333  )
334  {
335  Q.push(otherSpec);
336  this->activeSpecies_[otherSpec] = true;
337  }
338 
339  }
340  // second generation link only (for those without first link)
341  for (label v=0; v<NbrABInit2nd[u]; v++)
342  {
343  label otherSpec = rABOtherSpec2nd(u, v);
344  scalar rAB = (PAB2nd(u, v)+CAB2nd(u, v))/Den;
345  // the link is stronger than the user-defined tolerance
346  if
347  (
348  rAB >= this->tolerance()
349  && !this->activeSpecies_[otherSpec]
350  )
351  {
352  Q.push(otherSpec);
353  this->activeSpecies_[otherSpec] = true;
354  }
355  }
356  }
357  }
358 
360 }
361 
362 
363 // ************************************************************************* //
A FIFO stack based on a singly-linked list.
Definition: FIFOStack.H:51
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
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
scalar omega(const scalar p, const scalar T, const scalarField &c, const label li, scalar &omegaf, scalar &omegar) const
Net reaction rate.
Definition: Reaction.C:319
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const label nSpecie
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: reactionI.H:36
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:55
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, List< label > &ctos, DynamicList< label > &stoc, const label li)
Reduce the mechanism.
Definition: PFA.C:59
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))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
PFA(const IOdictionary &dict, chemistryModel< ThermoType > &chemistry)
Construct from components.
Definition: PFA.C:32
basicChemistryModel & chemistry
const multiComponentMixture< ThermoType > & mixture() const
Return reference to the mixture.
virtual ~PFA()
Destructor.
Definition: PFA.C:51
const volScalarField & T
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
Definition: reactionI.H:42
#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
An abstract class for methods of chemical mechanism reduction.
volScalarField & p
bool found