All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-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 "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_ = this->coeffsDict_.template lookup<scalar>("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 scalar p,
98  const scalar T,
99  const scalarField& c,
100  const label li
101 )
102 {
103  scalarField& completeC(this->chemistry_.completeC());
104  scalarField c1(this->chemistry_.nEqns(), 0.0);
105 
106  for (label i=0; i<this->nSpecie_; i++)
107  {
108  c1[i] = c[i];
109  completeC[i] = c[i];
110  }
111 
112  c1[this->nSpecie_] = T;
113  c1[this->nSpecie_+1] = p;
114 
115 
116  // Number of initialized rAB for each lines
117  Field<label> NbrABInit(this->nSpecie_,0);
118 
119  // Position of the initialized rAB, -1 when not initialized
120  RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
121  RectangularMatrix<scalar> CFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
122  RectangularMatrix<scalar> HFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
123  RectangularMatrix<scalar> OFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
124  RectangularMatrix<scalar> NFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
125  scalar CFlux(0.0), HFlux(0.0), OFlux(0.0), NFlux(0.0);
126  label nbPairs(0);
127 
128  // Index of the other species involved in the rABNum
129  RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
130 
131  scalar pf, cf, pr, cr;
132  label lRef, rRef;
133  forAll(this->chemistry_.reactions(), i)
134  {
135  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
136 
137  // for each reaction compute omegai
138  this->chemistry_.omega
139  (
140  R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef
141  );
142  scalar fr = mag(pf*cf)+mag(pr*cr);
143  scalar NCi(0.0),NHi(0.0),NOi(0.0),NNi(0.0);
144  forAll(R.lhs(),s)
145  {
146  label curIndex = R.lhs()[s].index;
147  scalar stoicCoeff = R.lhs()[s].stoichCoeff;
148  NCi += sC_[curIndex]*stoicCoeff;
149  NHi += sH_[curIndex]*stoicCoeff;
150  NOi += sO_[curIndex]*stoicCoeff;
151  NNi += sN_[curIndex]*stoicCoeff;
152  }
153  // element conservation means that total number of element is
154  // twice the number on the left
155  NCi *=2;
156  NHi *=2;
157  NOi *=2;
158  NNi *=2;
159  // then for each source/sink pairs, compute the element flux
160  forAll(R.lhs(), Ai)// compute the element flux
161  {
162  label A = R.lhs()[Ai].index;
163  scalar Acoeff = R.lhs()[Ai].stoichCoeff;
164 
165  forAll(R.rhs(), Bi)
166  {
167  label B = R.rhs()[Bi].index;
168  scalar Bcoeff = R.rhs()[Bi].stoichCoeff;
169  // if a pair in the reversed way has not been initialized
170  if (rABPos(B, A)==-1)
171  {
172  label otherS = rABPos(A, B);
173  if (otherS==-1)
174  {
175  rABPos(A, B) = NbrABInit[A]++;
176  otherS = rABPos(A, B);
177  nbPairs++;
178  rABOtherSpec(A, otherS) = B;
179  }
180  if (NCi>vSmall)
181  {
182  CFluxAB(A, otherS) +=
183  fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
184  }
185  if (NHi>vSmall)
186  {
187  HFluxAB(A, otherS) +=
188  fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
189  }
190  if (NOi>vSmall)
191  {
192  OFluxAB(A, otherS) +=
193  fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
194  }
195  if (NNi>vSmall)
196  {
197  NFluxAB(A, otherS) +=
198  fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
199  }
200  }
201  // If a pair BA is initialized,
202  // add the element flux to this pair
203  else
204  {
205  label otherS = rABPos(B, A);
206  if (NCi>vSmall)
207  {
208  CFluxAB(B, otherS) +=
209  fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
210  }
211  if (NHi>vSmall)
212  {
213  HFluxAB(B, otherS) +=
214  fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
215  }
216  if (NOi>vSmall)
217  {
218  OFluxAB(B, otherS) +=
219  fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
220  }
221  if (NNi>vSmall)
222  {
223  NFluxAB(B, otherS) +=
224  fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
225  }
226  }
227  if (NCi>vSmall)
228  {
229  CFlux += fr*Acoeff*sC_[A]*Bcoeff*sC_[B]/NCi;
230  }
231  if (NHi>vSmall)
232  {
233  HFlux += fr*Acoeff*sH_[A]*Bcoeff*sH_[B]/NHi;
234  }
235  if (NOi>vSmall)
236  {
237  OFlux += fr*Acoeff*sO_[A]*Bcoeff*sO_[B]/NOi;
238  }
239  if (NNi>vSmall)
240  {
241  NFlux += fr*Acoeff*sN_[A]*Bcoeff*sN_[B]/NNi;
242  }
243  }
244  }
245  }
246 
247  // Select species according to the total flux cutoff (1-tolerance)
248  // of the flux is included
249  label speciesNumber = 0;
250  for (label i=0; i<this->nSpecie_; i++)
251  {
252  this->activeSpecies_[i] = false;
253  }
254 
255  if (CFlux > vSmall)
256  {
257  SortableListEFA<scalar> pairsFlux(nbPairs);
258  labelList source(nbPairs);
259  labelList sink(nbPairs);
260  label nP(0);
261  for (int A=0; A<this->nSpecie_ ; A++)
262  {
263  for (int j=0; j<NbrABInit[A]; j++)
264  {
265  label pairIndex = nP++;
266  pairsFlux[pairIndex] = 0.0;
267  label B = rABOtherSpec(A, j);
268  pairsFlux[pairIndex] += CFluxAB(A, j);
269  source[pairIndex] = A;
270  sink[pairIndex] = B;
271  }
272  }
273 
274  // Sort in descending orders the source/sink pairs until the cutoff is
275  // reached. The sorting is done on part of the pairs because a small
276  // part of these pairs are responsible for a large part of the element
277  // flux
278  scalar cumFlux(0.0);
279  scalar threshold((1-this->tolerance())*CFlux);
280  label startPoint(0);
281  label nbToSort(static_cast<label> (nbPairs*sortPart_));
282  nbToSort = max(nbToSort,1);
283 
284  bool cumRespected(false);
285  while(!cumRespected)
286  {
287  if (startPoint >= nbPairs)
288  {
290  << "startPoint outside number of pairs without reaching"
291  << "100% flux"
292  << exit(FatalError);
293  }
294 
295  label nbi = min(nbToSort, nbPairs-startPoint);
296  pairsFlux.partialSort(nbi, startPoint);
297  const labelList& idx = pairsFlux.indices();
298  for (int i=0; i<nbi; i++)
299  {
300  cumFlux += pairsFlux[idx[startPoint+i]];
301  if (!this->activeSpecies_[source[idx[startPoint+i]]])
302  {
303  this->activeSpecies_[source[idx[startPoint+i]]] = true;
304  speciesNumber++;
305  }
306  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
307  {
308  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
309  speciesNumber++;
310  }
311  if (cumFlux >= threshold)
312  {
313  cumRespected = true;
314  break;
315  }
316  }
317  startPoint += nbToSort;
318  }
319  }
320 
321  if (HFlux > vSmall)
322  {
323  SortableListEFA<scalar> pairsFlux(nbPairs);
324  labelList source(nbPairs);
325  labelList sink(nbPairs);
326  label nP(0);
327  for (int A=0; A<this->nSpecie_ ; A++)
328  {
329  for (int j=0; j<NbrABInit[A]; j++)
330  {
331  label pairIndex = nP++;
332  pairsFlux[pairIndex] = 0.0;
333  label B = rABOtherSpec(A, j);
334  pairsFlux[pairIndex] += HFluxAB(A, j);
335  source[pairIndex] = A;
336  sink[pairIndex] = B;
337  }
338  }
339 
340  scalar cumFlux(0.0);
341  scalar threshold((1-this->tolerance())*HFlux);
342  label startPoint(0);
343  label nbToSort(static_cast<label> (nbPairs*sortPart_));
344  nbToSort = max(nbToSort,1);
345 
346  bool cumRespected(false);
347  while(!cumRespected)
348  {
349  if (startPoint >= nbPairs)
350  {
352  << "startPoint outside number of pairs without reaching"
353  << "100% flux"
354  << exit(FatalError);
355  }
356 
357  label nbi = min(nbToSort, nbPairs-startPoint);
358  pairsFlux.partialSort(nbi, startPoint);
359  const labelList& idx = pairsFlux.indices();
360  for (int i=0; i<nbi; i++)
361  {
362  cumFlux += pairsFlux[idx[startPoint+i]];
363 
364  if (!this->activeSpecies_[source[idx[startPoint+i]]])
365  {
366  this->activeSpecies_[source[idx[startPoint+i]]] = true;
367  speciesNumber++;
368  }
369  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
370  {
371  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
372  speciesNumber++;
373  }
374  if (cumFlux >= threshold)
375  {
376  cumRespected = true;
377  break;
378  }
379  }
380  startPoint += nbToSort;
381  }
382  }
383 
384  if (OFlux > vSmall)
385  {
386  SortableListEFA<scalar> pairsFlux(nbPairs);
387  labelList source(nbPairs);
388  labelList sink(nbPairs);
389  label nP(0);
390  for (int A=0; A<this->nSpecie_ ; A++)
391  {
392  for (int j=0; j<NbrABInit[A]; j++)
393  {
394  label pairIndex = nP++;
395  pairsFlux[pairIndex] = 0.0;
396  label B = rABOtherSpec(A, j);
397  pairsFlux[pairIndex] += OFluxAB(A, j);
398  source[pairIndex] = A;
399  sink[pairIndex] = B;
400  }
401  }
402  scalar cumFlux(0.0);
403  scalar threshold((1-this->tolerance())*OFlux);
404  label startPoint(0);
405  label nbToSort(static_cast<label> (nbPairs*sortPart_));
406  nbToSort = max(nbToSort,1);
407 
408  bool cumRespected(false);
409  while(!cumRespected)
410  {
411  if (startPoint >= nbPairs)
412  {
414  << "startPoint outside number of pairs without reaching"
415  << "100% flux"
416  << exit(FatalError);
417  }
418 
419  label nbi = min(nbToSort, nbPairs-startPoint);
420  pairsFlux.partialSort(nbi, startPoint);
421  const labelList& idx = pairsFlux.indices();
422  for (int i=0; i<nbi; i++)
423  {
424  cumFlux += pairsFlux[idx[startPoint+i]];
425 
426  if (!this->activeSpecies_[source[idx[startPoint+i]]])
427  {
428  this->activeSpecies_[source[idx[startPoint+i]]] = true;
429  speciesNumber++;
430  }
431  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
432  {
433  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
434  speciesNumber++;
435  }
436  if (cumFlux >= threshold)
437  {
438  cumRespected = true;
439  break;
440  }
441  }
442  startPoint += nbToSort;
443  }
444  }
445 
446  if (NFlux > vSmall)
447  {
448  SortableListEFA<scalar> pairsFlux(nbPairs);
449  labelList source(nbPairs);
450  labelList sink(nbPairs);
451  label nP(0);
452  for (int A=0; A<this->nSpecie_ ; A++)
453  {
454  for (int j=0; j<NbrABInit[A]; j++)
455  {
456  label pairIndex = nP++;
457  pairsFlux[pairIndex] = 0.0;
458  label B = rABOtherSpec(A, j);
459  pairsFlux[pairIndex] += NFluxAB(A, j);
460  source[pairIndex] = A;
461  sink[pairIndex] = B;
462  }
463  }
464  scalar cumFlux(0.0);
465  scalar threshold((1-this->tolerance())*NFlux);
466  label startPoint(0);
467  label nbToSort(static_cast<label> (nbPairs*sortPart_));
468  nbToSort = max(nbToSort,1);
469 
470  bool cumRespected(false);
471  while(!cumRespected)
472  {
473  if (startPoint >= nbPairs)
474  {
476  << "startPoint outside number of pairs without reaching"
477  << "100% flux"
478  << exit(FatalError);
479  }
480 
481  label nbi = min(nbToSort, nbPairs-startPoint);
482  pairsFlux.partialSort(nbi, startPoint);
483  const labelList& idx = pairsFlux.indices();
484  for (int i=0; i<nbi; i++)
485  {
486  cumFlux += pairsFlux[idx[startPoint+i]];
487 
488  if (!this->activeSpecies_[source[idx[startPoint+i]]])
489  {
490  this->activeSpecies_[source[idx[startPoint+i]]] = true;
491  speciesNumber++;
492  }
493  if (!this->activeSpecies_[sink[idx[startPoint+i]]])
494  {
495  this->activeSpecies_[sink[idx[startPoint+i]]] = true;
496  speciesNumber++;
497  }
498  if (cumFlux >= threshold)
499  {
500  cumRespected = true;
501  break;
502  }
503  }
504  startPoint += nbToSort;
505  }
506  }
507 
508  // Put a flag on the reactions containing at least one removed species
509  forAll(this->chemistry_.reactions(), i)
510  {
511  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
512  this->chemistry_.reactionsDisabled()[i] = false;
513  forAll(R.lhs(), s)
514  {
515  label ss = R.lhs()[s].index;
516  if (!this->activeSpecies_[ss])
517  {
518  this->chemistry_.reactionsDisabled()[i] = true;
519  break;
520  }
521  }
522  if (!this->chemistry_.reactionsDisabled()[i])
523  {
524  forAll(R.rhs(), s)
525  {
526  label ss = R.rhs()[s].index;
527  if (!this->activeSpecies_[ss])
528  {
529  this->chemistry_.reactionsDisabled()[i] = true;
530  break;
531  }
532  }
533  }
534  }
535 
536  this->NsSimp_ = speciesNumber;
537  scalarField& simplifiedC(this->chemistry_.simplifiedC());
538  simplifiedC.setSize(this->NsSimp_+2);
539  DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
540  s2c.setSize(this->NsSimp_);
541  Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
542  label j = 0;
543 
544  for (label i=0; i<this->nSpecie_; i++)
545  {
546  if (this->activeSpecies_[i])
547  {
548  s2c[j] = i;
549  simplifiedC[j] = c[i];
550  c2s[i] = j++;
551  if (!this->chemistry_.active(i))
552  {
553  this->chemistry_.setActive(i);
554  }
555  }
556  else
557  {
558  c2s[i] = -1;
559  }
560  }
561  simplifiedC[this->NsSimp_] = T;
562  simplifiedC[this->NsSimp_+1] = p;
563  this->chemistry_.setNsDAC(this->NsSimp_);
564  // change temporary Ns in chemistryModel
565  // to make the function nEqns working
566  this->chemistry_.setNSpecie(this->NsSimp_);
567 }
568 
569 
570 // ************************************************************************* //
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: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
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
const dimensionedScalar & c1
First radiation constant: default SI units: [W/m^2].
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: ReactionI.H:57
A list that is sorted upon construction or when explicitly requested with the sort() method...
BasicChemistryModel< rhoReactionThermo > & chemistry
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:175
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, const label li)
Reduce the mechanism.
Definition: EFA.C:96
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)
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
void partialSort(int M, int start)
Partial sort the list (if changed after construction time)
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