EFA.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 "EFA.H"
27 #include "SortableListEFA.H"
28 
29 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
30 
31 template<class ThermoType>
32 Foam::chemistryReductionMethods::EFA<ThermoType>::EFA
33 (
34  const IOdictionary& dict,
35  chemistryModel<ThermoType>& chemistry
36 )
37 :
38  chemistryReductionMethod<ThermoType>(dict, chemistry),
39  sC_(this->nSpecie(),0),
40  sH_(this->nSpecie(),0),
41  sO_(this->nSpecie(),0),
42  sN_(this->nSpecie(),0),
43  sortPart_(0.05)
44 {
45  for (label i=0; i<this->nSpecie(); i++)
46  {
47  const List<specieElement>& curSpecieComposition =
48  chemistry.mixture().specieComposition(i);
49 
50  // for all elements in the current species
51  forAll(curSpecieComposition, j)
52  {
53  const specieElement& curElement =
54  curSpecieComposition[j];
55  if (curElement.name() == "C")
56  {
57  sC_[i] = curElement.nAtoms();
58  }
59  else if (curElement.name() == "H")
60  {
61  sH_[i] = curElement.nAtoms();
62  }
63  else if (curElement.name() == "O")
64  {
65  sO_[i] = curElement.nAtoms();
66  }
67  else if (curElement.name() == "N")
68  {
69  sN_[i] = curElement.nAtoms();
70  }
71  else
72  {
73  Info<< "element not considered"<< endl;
74  }
75  }
76  }
77  if (this->coeffsDict_.found("sortPart"))
78  {
79  sortPart_ = this->coeffsDict_.template lookup<scalar>("sortPart");
80  }
81 }
82 
83 
84 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
85 
86 template<class ThermoType>
87 Foam::chemistryReductionMethods::EFA<ThermoType>::~EFA()
88 {}
89 
90 
91 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
92 
93 template<class ThermoType>
94 void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
95 (
96  const scalar p,
97  const scalar T,
98  const scalarField& c,
99  List<label>& ctos,
100  DynamicList<label>& stoc,
101  const label li
102 )
103 {
104  chemistryReductionMethod<ThermoType>::initReduceMechanism();
105 
106  scalarField c1(this->chemistry_.nEqns(), 0.0);
107 
108  for (label i=0; i<this->nSpecie(); i++)
109  {
110  c1[i] = c[i];
111  }
112 
113  c1[this->nSpecie()] = T;
114  c1[this->nSpecie()+1] = p;
115 
116 
117  // Number of initialised rAB for each lines
118  Field<label> NbrABInit(this->nSpecie(),0);
119 
120  // Position of the initialised rAB, -1 when not initialised
121  RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
122  RectangularMatrix<scalar> CFluxAB(this->nSpecie(), this->nSpecie(), 0.0);
123  RectangularMatrix<scalar> HFluxAB(this->nSpecie(), this->nSpecie(), 0.0);
124  RectangularMatrix<scalar> OFluxAB(this->nSpecie(), this->nSpecie(), 0.0);
125  RectangularMatrix<scalar> NFluxAB(this->nSpecie(), this->nSpecie(), 0.0);
126  scalar CFlux(0.0), HFlux(0.0), OFlux(0.0), NFlux(0.0);
127  label nbPairs(0);
128 
129  // Index of the other species involved in the rABNum
130  RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
131 
132  forAll(this->chemistry_.reactions(), i)
133  {
134  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
135 
136  // for each reaction compute omegai
137  scalar omegaf, omegar;
138  R.omega(p, T, c1, li, omegaf, omegar);
139 
140  scalar fr = mag(omegaf) + mag(omegar);
141  scalar NCi(0.0),NHi(0.0),NOi(0.0),NNi(0.0);
142  forAll(R.lhs(),s)
143  {
144  label curIndex = R.lhs()[s].index;
145  scalar stoicCoeff = R.lhs()[s].stoichCoeff;
146  NCi += sC_[curIndex]*stoicCoeff;
147  NHi += sH_[curIndex]*stoicCoeff;
148  NOi += sO_[curIndex]*stoicCoeff;
149  NNi += sN_[curIndex]*stoicCoeff;
150  }
151  // element conservation means that total number of element is
152  // twice the number on the left
153  NCi *=2;
154  NHi *=2;
155  NOi *=2;
156  NNi *=2;
157  // then for each source/sink pairs, compute the element flux
158  forAll(R.lhs(), Ai)// compute the element flux
159  {
160  label A = R.lhs()[Ai].index;
161  scalar Acoeff = R.lhs()[Ai].stoichCoeff;
162 
163  forAll(R.rhs(), Bi)
164  {
165  label B = R.rhs()[Bi].index;
166  scalar Bcoeff = R.rhs()[Bi].stoichCoeff;
167  // if a pair in the reversed way has not been initialised
168  if (rABPos(B, A)==-1)
169  {
170  label otherS = rABPos(A, B);
171  if (otherS==-1)
172  {
173  rABPos(A, B) = NbrABInit[A]++;
174  otherS = rABPos(A, B);
175  nbPairs++;
176  rABOtherSpec(A, otherS) = B;
177  }
178  if (NCi>vSmall)
179  {
180  CFluxAB(A, otherS) +=
181  fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
182  }
183  if (NHi>vSmall)
184  {
185  HFluxAB(A, otherS) +=
186  fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
187  }
188  if (NOi>vSmall)
189  {
190  OFluxAB(A, otherS) +=
191  fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
192  }
193  if (NNi>vSmall)
194  {
195  NFluxAB(A, otherS) +=
196  fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
197  }
198  }
199  // If a pair BA is initialised,
200  // add the element flux to this pair
201  else
202  {
203  label otherS = rABPos(B, A);
204  if (NCi>vSmall)
205  {
206  CFluxAB(B, otherS) +=
207  fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
208  }
209  if (NHi>vSmall)
210  {
211  HFluxAB(B, otherS) +=
212  fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
213  }
214  if (NOi>vSmall)
215  {
216  OFluxAB(B, otherS) +=
217  fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
218  }
219  if (NNi>vSmall)
220  {
221  NFluxAB(B, otherS) +=
222  fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
223  }
224  }
225  if (NCi>vSmall)
226  {
227  CFlux += fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
228  }
229  if (NHi>vSmall)
230  {
231  HFlux += fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
232  }
233  if (NOi>vSmall)
234  {
235  OFlux += fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
236  }
237  if (NNi>vSmall)
238  {
239  NFlux += fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
240  }
241  }
242  }
243  }
244 
245  // Select species according to the total flux cutoff (1-tolerance)
246  // of the flux is included
247  for (label i=0; i<this->nSpecie(); i++)
248  {
249  this->activeSpecies_[i] = false;
250  }
251 
252  if (CFlux > vSmall)
253  {
254  SortableListEFA<scalar> pairsFlux(nbPairs);
255  labelList source(nbPairs);
256  labelList sink(nbPairs);
257  label nP(0);
258  for (int A=0; A<this->nSpecie() ; A++)
259  {
260  for (int j=0; j<NbrABInit[A]; j++)
261  {
262  label pairIndex = nP++;
263  pairsFlux[pairIndex] = 0.0;
264  label B = rABOtherSpec(A, j);
265  pairsFlux[pairIndex] += CFluxAB(A, j);
266  source[pairIndex] = A;
267  sink[pairIndex] = B;
268  }
269  }
270 
271  // Sort in descending orders the source/sink pairs until the cutoff is
272  // reached. The sorting is done on part of the pairs because a small
273  // part of these pairs are responsible for a large part of the element
274  // flux
275  scalar cumFlux(0.0);
276  scalar threshold((1-this->tolerance())*CFlux);
277  label startPoint(0);
278  label nbToSort(static_cast<label> (nbPairs*sortPart_));
279  nbToSort = max(nbToSort,1);
280 
281  bool cumRespected(false);
282  while(!cumRespected)
283  {
284  if (startPoint >= nbPairs)
285  {
286  FatalErrorInFunction
287  << "startPoint outside number of pairs without reaching"
288  << "100% flux"
289  << exit(FatalError);
290  }
291 
292  label nbi = min(nbToSort, nbPairs-startPoint);
293  pairsFlux.partialSort(nbi, startPoint);
294  const labelList& idx = pairsFlux.indices();
295  for (int i=0; i<nbi; i++)
296  {
297  cumFlux += pairsFlux[idx[startPoint+i]];
298  if (!this->activeSpecies_[source[idx[startPoint+i]]])
299  {
300  this->activeSpecies_[source[idx[startPoint+i]]] = true;
301  }
302  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
303  {
304  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
305  }
306  if (cumFlux >= threshold)
307  {
308  cumRespected = true;
309  break;
310  }
311  }
312  startPoint += nbToSort;
313  }
314  }
315 
316  if (HFlux > vSmall)
317  {
318  SortableListEFA<scalar> pairsFlux(nbPairs);
319  labelList source(nbPairs);
320  labelList sink(nbPairs);
321  label nP(0);
322  for (int A=0; A<this->nSpecie() ; A++)
323  {
324  for (int j=0; j<NbrABInit[A]; j++)
325  {
326  label pairIndex = nP++;
327  pairsFlux[pairIndex] = 0.0;
328  label B = rABOtherSpec(A, j);
329  pairsFlux[pairIndex] += HFluxAB(A, j);
330  source[pairIndex] = A;
331  sink[pairIndex] = B;
332  }
333  }
334 
335  scalar cumFlux(0.0);
336  scalar threshold((1-this->tolerance())*HFlux);
337  label startPoint(0);
338  label nbToSort(static_cast<label> (nbPairs*sortPart_));
339  nbToSort = max(nbToSort,1);
340 
341  bool cumRespected(false);
342  while(!cumRespected)
343  {
344  if (startPoint >= nbPairs)
345  {
346  FatalErrorInFunction
347  << "startPoint outside number of pairs without reaching"
348  << "100% flux"
349  << exit(FatalError);
350  }
351 
352  label nbi = min(nbToSort, nbPairs-startPoint);
353  pairsFlux.partialSort(nbi, startPoint);
354  const labelList& idx = pairsFlux.indices();
355  for (int i=0; i<nbi; i++)
356  {
357  cumFlux += pairsFlux[idx[startPoint+i]];
358 
359  if (!this->activeSpecies_[source[idx[startPoint+i]]])
360  {
361  this->activeSpecies_[source[idx[startPoint+i]]] = true;
362  }
363  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
364  {
365  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
366  }
367  if (cumFlux >= threshold)
368  {
369  cumRespected = true;
370  break;
371  }
372  }
373  startPoint += nbToSort;
374  }
375  }
376 
377  if (OFlux > vSmall)
378  {
379  SortableListEFA<scalar> pairsFlux(nbPairs);
380  labelList source(nbPairs);
381  labelList sink(nbPairs);
382  label nP(0);
383  for (int A=0; A<this->nSpecie() ; A++)
384  {
385  for (int j=0; j<NbrABInit[A]; j++)
386  {
387  label pairIndex = nP++;
388  pairsFlux[pairIndex] = 0.0;
389  label B = rABOtherSpec(A, j);
390  pairsFlux[pairIndex] += OFluxAB(A, j);
391  source[pairIndex] = A;
392  sink[pairIndex] = B;
393  }
394  }
395  scalar cumFlux(0.0);
396  scalar threshold((1-this->tolerance())*OFlux);
397  label startPoint(0);
398  label nbToSort(static_cast<label> (nbPairs*sortPart_));
399  nbToSort = max(nbToSort,1);
400 
401  bool cumRespected(false);
402  while(!cumRespected)
403  {
404  if (startPoint >= nbPairs)
405  {
406  FatalErrorInFunction
407  << "startPoint outside number of pairs without reaching"
408  << "100% flux"
409  << exit(FatalError);
410  }
411 
412  label nbi = min(nbToSort, nbPairs-startPoint);
413  pairsFlux.partialSort(nbi, startPoint);
414  const labelList& idx = pairsFlux.indices();
415  for (int i=0; i<nbi; i++)
416  {
417  cumFlux += pairsFlux[idx[startPoint+i]];
418 
419  if (!this->activeSpecies_[source[idx[startPoint+i]]])
420  {
421  this->activeSpecies_[source[idx[startPoint+i]]] = true;
422  }
423  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
424  {
425  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
426  }
427  if (cumFlux >= threshold)
428  {
429  cumRespected = true;
430  break;
431  }
432  }
433  startPoint += nbToSort;
434  }
435  }
436 
437  if (NFlux > vSmall)
438  {
439  SortableListEFA<scalar> pairsFlux(nbPairs);
440  labelList source(nbPairs);
441  labelList sink(nbPairs);
442  label nP(0);
443  for (int A=0; A<this->nSpecie() ; A++)
444  {
445  for (int j=0; j<NbrABInit[A]; j++)
446  {
447  label pairIndex = nP++;
448  pairsFlux[pairIndex] = 0.0;
449  label B = rABOtherSpec(A, j);
450  pairsFlux[pairIndex] += NFluxAB(A, j);
451  source[pairIndex] = A;
452  sink[pairIndex] = B;
453  }
454  }
455  scalar cumFlux(0.0);
456  scalar threshold((1-this->tolerance())*NFlux);
457  label startPoint(0);
458  label nbToSort(static_cast<label> (nbPairs*sortPart_));
459  nbToSort = max(nbToSort,1);
460 
461  bool cumRespected(false);
462  while(!cumRespected)
463  {
464  if (startPoint >= nbPairs)
465  {
466  FatalErrorInFunction
467  << "startPoint outside number of pairs without reaching"
468  << "100% flux"
469  << exit(FatalError);
470  }
471 
472  label nbi = min(nbToSort, nbPairs-startPoint);
473  pairsFlux.partialSort(nbi, startPoint);
474  const labelList& idx = pairsFlux.indices();
475  for (int i=0; i<nbi; i++)
476  {
477  cumFlux += pairsFlux[idx[startPoint+i]];
478 
479  if (!this->activeSpecies_[source[idx[startPoint+i]]])
480  {
481  this->activeSpecies_[source[idx[startPoint+i]]] = true;
482  }
483  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
484  {
485  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
486  }
487  if (cumFlux >= threshold)
488  {
489  cumRespected = true;
490  break;
491  }
492  }
493  startPoint += nbToSort;
494  }
495  }
496 
497  chemistryReductionMethod<ThermoType>::endReduceMechanism(ctos, stoc);
498 }
499 
500 
501 // ************************************************************************* //
A
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
B
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Foam::IOdictionary
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:91
Foam::Reaction
Simple extension of ThermoType to handle reaction kinetics in addition to the equilibrium thermodynam...
Definition: Reaction.H:72
Foam::RectangularMatrix
A templated 2D rectangular m x n matrix of objects of <Type>.
Definition: RectangularMatrix.H:56
Foam::SortableListEFA
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableListEFA.H:53
Foam::SortableListEFA::indices
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableListEFA.H:102
Foam::SortableListEFA::partialSort
void partialSort(int M, int start)
Partial sort the list (if changed after construction time)
Definition: SortableListEFA.C:112
Foam::chemistryModel
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
Definition: chemistryModel.H:90
Foam::chemistryReductionMethod
An abstract class for methods of chemical mechanism reduction.
Definition: chemistryReductionMethod.H:60
Foam::chemistryReductionMethod::initReduceMechanism
void initReduceMechanism()
Protected Member Functions.
Definition: chemistryReductionMethod.C:88
Foam::chemistryReductionMethod::coeffsDict_
const dictionary coeffsDict_
Dictionary that store the algorithm data.
Definition: chemistryReductionMethod.H:66
Foam::chemistryReductionMethod::nSpecie
label nSpecie()
Return the number of species.
Definition: chemistryReductionMethodI.H:31
Foam::chemistryReductionMethod::endReduceMechanism
void endReduceMechanism(List< label > &ctos, DynamicList< label > &stoc)
End reduction of the mechanism.
Definition: chemistryReductionMethod.C:99
Foam::chemistryReductionMethods::EFA::~EFA
virtual ~EFA()
Destructor.
Definition: EFA.C:87
Foam::chemistryReductionMethods::EFA::EFA
EFA(const IOdictionary &dict, chemistryModel< ThermoType > &chemistry)
Construct from components.
Definition: EFA.C:33
Foam::chemistryReductionMethods::EFA::reduceMechanism
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: EFA.C:95
Foam::dictionary::found
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:509
Foam::specieElement::nAtoms
label nAtoms() const
Return the number of atoms of this element in the specie.
Definition: specieElementI.H:62
Foam::specieElement::name
const word & name() const
Return the name of the element.
Definition: specieElementI.H:50
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
s
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.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Foam::constant::physicoChemical::c1
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Foam::label
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
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
Foam::Info
messageStream Info
Foam::min
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Definition: fvMeshStitchersMoving.C:110
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
Foam::max
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Definition: fvMeshStitchersMoving.C:104
Foam::R
static scalar R(const scalar a, const scalar x)
Definition: invIncGamma.C:102
Foam::FatalError
error FatalError
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:55
nSpecie
const label nSpecie
Definition: readInitialConditions.H:24
dict
dictionary dict
Definition: searchingEngine.H:14
p
volScalarField & p
Definition: createFieldRefs.H:4
chemistry
basicChemistryModel & chemistry
Definition: createFieldRefs.H:1