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-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 "EFA.H"
27 #include "SortableListEFA.H"
28 
29 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
30 
31 template<class CompType, class ThermoType>
33 (
34  const IOdictionary& dict,
36 )
37 :
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  const List<List<specieElement>>& specieComposition =
46  this->chemistry_.specieComp();
47  for (label i=0; i<this->nSpecie_; i++)
48  {
49  const List<specieElement>& curSpecieComposition =
50  specieComposition[i];
51  // for all elements in the current species
52  forAll(curSpecieComposition, j)
53  {
54  const specieElement& curElement =
55  curSpecieComposition[j];
56  if (curElement.name() == "C")
57  {
58  sC_[i] = curElement.nAtoms();
59  }
60  else if (curElement.name() == "H")
61  {
62  sH_[i] = curElement.nAtoms();
63  }
64  else if (curElement.name() == "O")
65  {
66  sO_[i] = curElement.nAtoms();
67  }
68  else if (curElement.name() == "N")
69  {
70  sN_[i] = curElement.nAtoms();
71  }
72  else
73  {
74  Info<< "element not considered"<< endl;
75  }
76  }
77  }
78  if (this->coeffsDict_.found("sortPart"))
79  {
80  sortPart_ = readScalar(this->coeffsDict_.lookup("sortPart"));
81  }
82 }
83 
84 
85 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
86 
87 template<class CompType, class ThermoType>
89 {}
90 
91 
92 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
93 
94 template<class CompType, class ThermoType>
96 (
97  const scalarField &c,
98  const scalar T,
99  const scalar p
100 )
101 {
102  scalarField& completeC(this->chemistry_.completeC());
103  scalarField c1(this->chemistry_.nEqns(), 0.0);
104 
105  for (label i=0; i<this->nSpecie_; i++)
106  {
107  c1[i] = c[i];
108  completeC[i] = c[i];
109  }
110 
111  c1[this->nSpecie_] = T;
112  c1[this->nSpecie_+1] = p;
113 
114 
115  // Number of initialized rAB for each lines
116  Field<label> NbrABInit(this->nSpecie_,0);
117 
118  // Position of the initialized rAB, -1 when not initialized
119  RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
120  RectangularMatrix<scalar> CFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
121  RectangularMatrix<scalar> HFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
122  RectangularMatrix<scalar> OFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
123  RectangularMatrix<scalar> NFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
124  scalar CFlux(0.0), HFlux(0.0), OFlux(0.0), NFlux(0.0);
125  label nbPairs(0);
126 
127  // Index of the other species involved in the rABNum
128  RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
129 
130  scalar pf, cf, pr, cr;
131  label lRef, rRef;
132  forAll(this->chemistry_.reactions(), i)
133  {
134  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
135  // for each reaction compute omegai
136  this->chemistry_.omega
137  (
138  R, c1, T, p, pf, cf, lRef, pr, cr, rRef
139  );
140  scalar fr = mag(pf*cf)+mag(pr*cr);
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 initialized
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 initialized,
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  label speciesNumber = 0;
248  for (label i=0; i<this->nSpecie_; i++)
249  {
250  this->activeSpecies_[i] = false;
251  }
252 
253  if (CFlux > vSmall)
254  {
255  SortableListEFA<scalar> pairsFlux(nbPairs);
256  labelList source(nbPairs);
257  labelList sink(nbPairs);
258  label nP(0);
259  for (int A=0; A<this->nSpecie_ ; A++)
260  {
261  for (int j=0; j<NbrABInit[A]; j++)
262  {
263  label pairIndex = nP++;
264  pairsFlux[pairIndex] = 0.0;
265  label B = rABOtherSpec(A, j);
266  pairsFlux[pairIndex] += CFluxAB(A, j);
267  source[pairIndex] = A;
268  sink[pairIndex] = B;
269  }
270  }
271 
272  // Sort in descending orders the source/sink pairs until the cutoff is
273  // reached. The sorting is done on part of the pairs because a small
274  // part of these pairs are responsible for a large part of the element
275  // flux
276  scalar cumFlux(0.0);
277  scalar threshold((1-this->tolerance())*CFlux);
278  label startPoint(0);
279  label nbToSort(static_cast<label> (nbPairs*sortPart_));
280  nbToSort = max(nbToSort,1);
281 
282  bool cumRespected(false);
283  while(!cumRespected)
284  {
285  if (startPoint >= nbPairs)
286  {
288  << "startPoint outside number of pairs without reaching"
289  << "100% flux"
290  << exit(FatalError);
291  }
292 
293  label nbi = min(nbToSort, nbPairs-startPoint);
294  pairsFlux.partialSort(nbi, startPoint);
295  const labelList& idx = pairsFlux.indices();
296  for (int i=0; i<nbi; i++)
297  {
298  cumFlux += pairsFlux[idx[startPoint+i]];
299  if (!this->activeSpecies_[source[idx[startPoint+i]]])
300  {
301  this->activeSpecies_[source[idx[startPoint+i]]] = true;
302  speciesNumber++;
303  }
304  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
305  {
306  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
307  speciesNumber++;
308  }
309  if (cumFlux >= threshold)
310  {
311  cumRespected = true;
312  break;
313  }
314  }
315  startPoint += nbToSort;
316  }
317  }
318 
319  if (HFlux > vSmall)
320  {
321  SortableListEFA<scalar> pairsFlux(nbPairs);
322  labelList source(nbPairs);
323  labelList sink(nbPairs);
324  label nP(0);
325  for (int A=0; A<this->nSpecie_ ; A++)
326  {
327  for (int j=0; j<NbrABInit[A]; j++)
328  {
329  label pairIndex = nP++;
330  pairsFlux[pairIndex] = 0.0;
331  label B = rABOtherSpec(A, j);
332  pairsFlux[pairIndex] += HFluxAB(A, j);
333  source[pairIndex] = A;
334  sink[pairIndex] = B;
335  }
336  }
337 
338  scalar cumFlux(0.0);
339  scalar threshold((1-this->tolerance())*HFlux);
340  label startPoint(0);
341  label nbToSort(static_cast<label> (nbPairs*sortPart_));
342  nbToSort = max(nbToSort,1);
343 
344  bool cumRespected(false);
345  while(!cumRespected)
346  {
347  if (startPoint >= nbPairs)
348  {
350  << "startPoint outside number of pairs without reaching"
351  << "100% flux"
352  << exit(FatalError);
353  }
354 
355  label nbi = min(nbToSort, nbPairs-startPoint);
356  pairsFlux.partialSort(nbi, startPoint);
357  const labelList& idx = pairsFlux.indices();
358  for (int i=0; i<nbi; i++)
359  {
360  cumFlux += pairsFlux[idx[startPoint+i]];
361 
362  if (!this->activeSpecies_[source[idx[startPoint+i]]])
363  {
364  this->activeSpecies_[source[idx[startPoint+i]]] = true;
365  speciesNumber++;
366  }
367  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
368  {
369  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
370  speciesNumber++;
371  }
372  if (cumFlux >= threshold)
373  {
374  cumRespected = true;
375  break;
376  }
377  }
378  startPoint += nbToSort;
379  }
380  }
381 
382  if (OFlux > vSmall)
383  {
384  SortableListEFA<scalar> pairsFlux(nbPairs);
385  labelList source(nbPairs);
386  labelList sink(nbPairs);
387  label nP(0);
388  for (int A=0; A<this->nSpecie_ ; A++)
389  {
390  for (int j=0; j<NbrABInit[A]; j++)
391  {
392  label pairIndex = nP++;
393  pairsFlux[pairIndex] = 0.0;
394  label B = rABOtherSpec(A, j);
395  pairsFlux[pairIndex] += OFluxAB(A, j);
396  source[pairIndex] = A;
397  sink[pairIndex] = B;
398  }
399  }
400  scalar cumFlux(0.0);
401  scalar threshold((1-this->tolerance())*OFlux);
402  label startPoint(0);
403  label nbToSort(static_cast<label> (nbPairs*sortPart_));
404  nbToSort = max(nbToSort,1);
405 
406  bool cumRespected(false);
407  while(!cumRespected)
408  {
409  if (startPoint >= nbPairs)
410  {
412  << "startPoint outside number of pairs without reaching"
413  << "100% flux"
414  << exit(FatalError);
415  }
416 
417  label nbi = min(nbToSort, nbPairs-startPoint);
418  pairsFlux.partialSort(nbi, startPoint);
419  const labelList& idx = pairsFlux.indices();
420  for (int i=0; i<nbi; i++)
421  {
422  cumFlux += pairsFlux[idx[startPoint+i]];
423 
424  if (!this->activeSpecies_[source[idx[startPoint+i]]])
425  {
426  this->activeSpecies_[source[idx[startPoint+i]]] = true;
427  speciesNumber++;
428  }
429  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
430  {
431  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
432  speciesNumber++;
433  }
434  if (cumFlux >= threshold)
435  {
436  cumRespected = true;
437  break;
438  }
439  }
440  startPoint += nbToSort;
441  }
442  }
443 
444  if (NFlux > vSmall)
445  {
446  SortableListEFA<scalar> pairsFlux(nbPairs);
447  labelList source(nbPairs);
448  labelList sink(nbPairs);
449  label nP(0);
450  for (int A=0; A<this->nSpecie_ ; A++)
451  {
452  for (int j=0; j<NbrABInit[A]; j++)
453  {
454  label pairIndex = nP++;
455  pairsFlux[pairIndex] = 0.0;
456  label B = rABOtherSpec(A, j);
457  pairsFlux[pairIndex] += NFluxAB(A, j);
458  source[pairIndex] = A;
459  sink[pairIndex] = B;
460  }
461  }
462  scalar cumFlux(0.0);
463  scalar threshold((1-this->tolerance())*NFlux);
464  label startPoint(0);
465  label nbToSort(static_cast<label> (nbPairs*sortPart_));
466  nbToSort = max(nbToSort,1);
467 
468  bool cumRespected(false);
469  while(!cumRespected)
470  {
471  if (startPoint >= nbPairs)
472  {
474  << "startPoint outside number of pairs without reaching"
475  << "100% flux"
476  << exit(FatalError);
477  }
478 
479  label nbi = min(nbToSort, nbPairs-startPoint);
480  pairsFlux.partialSort(nbi, startPoint);
481  const labelList& idx = pairsFlux.indices();
482  for (int i=0; i<nbi; i++)
483  {
484  cumFlux += pairsFlux[idx[startPoint+i]];
485 
486  if (!this->activeSpecies_[source[idx[startPoint+i]]])
487  {
488  this->activeSpecies_[source[idx[startPoint+i]]] = true;
489  speciesNumber++;
490  }
491  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
492  {
493  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
494  speciesNumber++;
495  }
496  if (cumFlux >= threshold)
497  {
498  cumRespected = true;
499  break;
500  }
501  }
502  startPoint += nbToSort;
503  }
504  }
505 
506  // Put a flag on the reactions containing at least one removed species
507  forAll(this->chemistry_.reactions(), i)
508  {
509  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
510  this->chemistry_.reactionsDisabled()[i] = false;
511  forAll(R.lhs(), s)
512  {
513  label ss = R.lhs()[s].index;
514  if (!this->activeSpecies_[ss])
515  {
516  this->chemistry_.reactionsDisabled()[i] = true;
517  break;
518  }
519  }
520  if (!this->chemistry_.reactionsDisabled()[i])
521  {
522  forAll(R.rhs(), s)
523  {
524  label ss = R.rhs()[s].index;
525  if (!this->activeSpecies_[ss])
526  {
527  this->chemistry_.reactionsDisabled()[i] = true;
528  break;
529  }
530  }
531  }
532  }
533 
534  this->NsSimp_ = speciesNumber;
535  scalarField& simplifiedC(this->chemistry_.simplifiedC());
536  simplifiedC.setSize(this->NsSimp_+2);
537  DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
538  s2c.setSize(this->NsSimp_);
539  Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
540  label j = 0;
541 
542  for (label i=0; i<this->nSpecie_; i++)
543  {
544  if (this->activeSpecies_[i])
545  {
546  s2c[j] = i;
547  simplifiedC[j] = c[i];
548  c2s[i] = j++;
549  if (!this->chemistry_.active(i))
550  {
551  this->chemistry_.setActive(i);
552  }
553  }
554  else
555  {
556  c2s[i] = -1;
557  }
558  }
559  simplifiedC[this->NsSimp_] = T;
560  simplifiedC[this->NsSimp_+1] = p;
561  this->chemistry_.setNsDAC(this->NsSimp_);
562  // change temporary Ns in chemistryModel
563  // to make the function nEqns working
564  this->chemistry_.setNSpecie(this->NsSimp_);
565 }
566 
567 
568 // ************************************************************************* //
#define readScalar
Definition: doubleScalar.C:38
EFA(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
Definition: EFA.C:33
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: ReactionI.H:58
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:60
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
A list that is sorted upon construction or when explicitly requested with the sort() method...
BasicChemistryModel< rhoReactionThermo > & chemistry
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:55
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
Definition: EFA.C:96
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:52
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:163
const word & name() const
Return the name of the element.
label nAtoms() const
Return the number of atoms of this element in the specie.
A templated 2D rectangular m x n matrix of objects of <Type>.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const volScalarField & T
void setSize(const label)
Reset size of List.
Definition: List.C:281
virtual ~EFA()
Destructor.
Definition: EFA.C:88
#define R(A, B, C, D, E, F, K, M)
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
void partialSort(int M, int start)
Partial sort the list (if changed after construction time)
An abstract class for methods of chemical mechanism reduction.
volScalarField & p
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
Definition: ReactionI.H:66