EFA.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "EFA.H"
29#include "SortableListEFA.H"
30
31// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32
33template<class CompType, class ThermoType>
34Foam::chemistryReductionMethods::EFA<CompType, ThermoType>::EFA
35(
36 const IOdictionary& dict,
37 TDACChemistryModel<CompType, ThermoType>& chemistry
38)
39:
40 chemistryReductionMethod<CompType, ThermoType>(dict, chemistry),
41 sC_(this->nSpecie_, Zero),
42 sH_(this->nSpecie_, Zero),
43 sO_(this->nSpecie_, Zero),
44 sN_(this->nSpecie_, Zero),
45 sortPart_(0.05)
46{
47 const List<List<specieElement>>& specieComposition =
48 this->chemistry_.specieComp();
49 for (label i=0; i<this->nSpecie_; i++)
50 {
51 const List<specieElement>& curSpecieComposition =
52 specieComposition[i];
53 // for all elements in the current species
54 forAll(curSpecieComposition, j)
55 {
56 const specieElement& curElement =
57 curSpecieComposition[j];
58 if (curElement.name() == "C")
59 {
60 sC_[i] = curElement.nAtoms();
61 }
62 else if (curElement.name() == "H")
63 {
64 sH_[i] = curElement.nAtoms();
65 }
66 else if (curElement.name() == "O")
67 {
68 sO_[i] = curElement.nAtoms();
69 }
70 else if (curElement.name() == "N")
71 {
72 sN_[i] = curElement.nAtoms();
73 }
74 else
75 {
76 Info<< "element not considered"<< endl;
77 }
78 }
79 }
80
81 this->coeffsDict_.readIfPresent("sortPart", sortPart_);
82}
83
84
85// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
86
87template<class CompType, class ThermoType>
88Foam::chemistryReductionMethods::EFA<CompType, ThermoType>::~EFA()
89{}
90
91
92// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
93
94template<class CompType, class ThermoType>
95void Foam::chemistryReductionMethods::EFA<CompType, ThermoType>::reduceMechanism
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(), Zero);
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_, Zero);
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_, Zero);
121 RectangularMatrix<scalar> HFluxAB(this->nSpecie_, this->nSpecie_, Zero);
122 RectangularMatrix<scalar> OFluxAB(this->nSpecie_, this->nSpecie_, Zero);
123 RectangularMatrix<scalar> NFluxAB(this->nSpecie_, this->nSpecie_, Zero);
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 {
287 FatalErrorInFunction
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 {
349 FatalErrorInFunction
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 {
411 FatalErrorInFunction
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 {
473 FatalErrorInFunction
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// ************************************************************************* //
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
B
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
R
#define R(A, B, C, D, E, F, K, M)
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
Foam::DynamicList::setSize
void setSize(const label n)
Same as resize()
Definition: DynamicList.H:221
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:77
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:218
Foam::Reaction
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:73
Foam::RectangularMatrix
A templated (M x N) rectangular matrix of objects of <Type>, containing M*N elements,...
Definition: RectangularMatrix.H:60
Foam::SortableListEFA
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableListEFA.H:55
Foam::SortableListEFA::indices
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableListEFA.H:104
Foam::SortableListEFA::partialSort
void partialSort(int M, int start)
Partial sort the list (if changed after construction time)
Definition: SortableListEFA.C:114
Foam::TDACChemistryModel
Extends StandardChemistryModel by adding the TDAC method.
Definition: TDACChemistryModel.H:86
Foam::TDACChemistryModel::specieComp
List< List< specieElement > > & specieComp()
Definition: TDACChemistryModelI.H:168
Foam::chemistryReductionMethod
An abstract class for methods of chemical mechanism reduction.
Definition: chemistryReductionMethod.H:57
Foam::chemistryReductionMethod::chemistry_
TDACChemistryModel< CompType, ThermoType > & chemistry_
Definition: chemistryReductionMethod.H:72
Foam::chemistryReductionMethod::nSpecie_
const label nSpecie_
Number of species.
Definition: chemistryReductionMethod.H:81
Foam::chemistryReductionMethod::coeffsDict_
const dictionary coeffsDict_
Dictionary that store the algorithm data.
Definition: chemistryReductionMethod.H:64
Foam::chemistryReductionMethods::EFA::~EFA
virtual ~EFA()
Destructor.
Definition: EFA.C:88
Foam::chemistryReductionMethods::EFA::reduceMechanism
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
Definition: EFA.C:96
Foam::dictionary::readIfPresent
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:405
Foam::specieElement::nAtoms
label nAtoms() const
Return the number of atoms of this element in the specie.
Definition: specieElementI.H:68
Foam::specieElement::name
const word & name() const
Return the name of the element.
Definition: specieElementI.H:56
p
volScalarField & p
Definition: createFieldRefs.H:8
chemistry
BasicChemistryModel< psiReactionThermo > & chemistry
Definition: createFieldRefs.H:1
T
const volScalarField & T
Definition: createFieldRefs.H:2
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::FatalError
error FatalError
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dict
dictionary dict
Definition: searchingEngine.H:14
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333