TDACChemistryModel.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-2017 OpenFOAM Foundation
9  Copyright (C) 2016-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "TDACChemistryModel.H"
30 #include "UniformField.H"
31 #include "localEulerDdtScheme.H"
32 #include "clockTime.H"
33 
34 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
35 
36 template<class ReactionThermo, class ThermoType>
38 (
39  ReactionThermo& thermo
40 )
41 :
43  variableTimeStep_
44  (
45  this->mesh().time().controlDict().getOrDefault
46  (
47  "adjustTimeStep",
48  false
49  )
50  || fv::localEulerDdt::enabled(this->mesh())
51  ),
52  timeSteps_(0),
53  NsDAC_(this->nSpecie_),
54  completeC_(this->nSpecie_, 0),
55  reactionsDisabled_(this->reactions_.size(), false),
56  specieComp_(this->nSpecie_),
57  completeToSimplifiedIndex_(this->nSpecie_, -1),
58  simplifiedToCompleteIndex_(this->nSpecie_),
59  tabulationResults_
60  (
61  IOobject
62  (
63  thermo.phasePropertyName("TabulationResults"),
64  this->time().timeName(),
65  this->mesh(),
66  IOobject::NO_READ,
67  IOobject::AUTO_WRITE
68  ),
69  this->mesh(),
71  )
72 {
73  basicSpecieMixture& composition = this->thermo().composition();
74 
75  // Store the species composition according to the species index
76  speciesTable speciesTab = composition.species();
77 
79  (
81  .specieComposition()
82  );
83 
84  forAll(specieComp_, i)
85  {
86  specieComp_[i] = (specCompPtr.ref())[this->Y()[i].member()];
87  }
88 
90  (
91  *this,
92  *this
93  );
94 
95  // When the mechanism reduction method is used, the 'active' flag for every
96  // species should be initialized (by default 'active' is true)
97  if (mechRed_->active())
98  {
99  forAll(this->Y(), i)
100  {
101  IOobject header
102  (
103  this->Y()[i].name(),
104  this->mesh().time().timeName(),
105  this->mesh(),
106  IOobject::NO_READ
107  );
108 
109  // Check if the species file is provided, if not set inactive
110  // and NO_WRITE
111  if (!header.typeHeaderOk<volScalarField>(true))
112  {
113  composition.setInactive(i);
114  this->Y()[i].writeOpt(IOobject::NO_WRITE);
115  }
116  }
117  }
118 
120  (
121  *this,
122  *this
123  );
124 
125  if (mechRed_->log())
126  {
127  cpuReduceFile_ = logFile("cpu_reduce.out");
128  nActiveSpeciesFile_ = logFile("nActiveSpecies.out");
129  }
130 
131  if (tabulation_->log())
132  {
133  cpuAddFile_ = logFile("cpu_add.out");
134  cpuGrowFile_ = logFile("cpu_grow.out");
135  cpuRetrieveFile_ = logFile("cpu_retrieve.out");
136  }
137 
138  if (mechRed_->log() || tabulation_->log())
139  {
140  cpuSolveFile_ = logFile("cpu_solve.out");
141  }
142 }
143 
144 
145 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
146 
147 template<class ReactionThermo, class ThermoType>
149 {}
150 
151 
152 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
153 
154 template<class ReactionThermo, class ThermoType>
156 (
157  const scalarField& c, // Contains all species even when mechRed is active
158  const scalar T,
159  const scalar p,
160  scalarField& dcdt
161 ) const
162 {
163  const bool reduced = mechRed_->active();
164 
165  scalar pf, cf, pr, cr;
166  label lRef, rRef;
167 
168  dcdt = Zero;
169 
170  forAll(this->reactions_, i)
171  {
172  if (!reactionsDisabled_[i])
173  {
174  const Reaction<ThermoType>& R = this->reactions_[i];
175 
176  scalar omegai = omega
177  (
178  R, c, T, p, pf, cf, lRef, pr, cr, rRef
179  );
180 
181  forAll(R.lhs(), s)
182  {
183  label si = R.lhs()[s].index;
184  if (reduced)
185  {
186  si = completeToSimplifiedIndex_[si];
187  }
188 
189  const scalar sl = R.lhs()[s].stoichCoeff;
190  dcdt[si] -= sl*omegai;
191  }
192  forAll(R.rhs(), s)
193  {
194  label si = R.rhs()[s].index;
195  if (reduced)
196  {
197  si = completeToSimplifiedIndex_[si];
198  }
199 
200  const scalar sr = R.rhs()[s].stoichCoeff;
201  dcdt[si] += sr*omegai;
202  }
203  }
204  }
205 }
206 
207 
208 template<class ReactionThermo, class ThermoType>
210 (
211  const Reaction<ThermoType>& R,
212  const scalarField& c, // Contains all species even when mechRed is active
213  const scalar T,
214  const scalar p,
215  scalar& pf,
216  scalar& cf,
217  label& lRef,
218  scalar& pr,
219  scalar& cr,
220  label& rRef
221 ) const
222 {
223  const scalar kf = R.kf(p, T, c);
224  const scalar kr = R.kr(kf, p, T, c);
225 
226  const label Nl = R.lhs().size();
227  const label Nr = R.rhs().size();
228 
229  label slRef = 0;
230  lRef = R.lhs()[slRef].index;
231 
232  pf = kf;
233  for (label s=1; s<Nl; s++)
234  {
235  const label si = R.lhs()[s].index;
236 
237  if (c[si] < c[lRef])
238  {
239  const scalar exp = R.lhs()[slRef].exponent;
240  pf *= pow(max(c[lRef], 0.0), exp);
241  lRef = si;
242  slRef = s;
243  }
244  else
245  {
246  const scalar exp = R.lhs()[s].exponent;
247  pf *= pow(max(c[si], 0.0), exp);
248  }
249  }
250  cf = max(c[lRef], 0.0);
251 
252  {
253  const scalar exp = R.lhs()[slRef].exponent;
254  if (exp < 1)
255  {
256  if (cf > SMALL)
257  {
258  pf *= pow(cf, exp - 1);
259  }
260  else
261  {
262  pf = 0;
263  }
264  }
265  else
266  {
267  pf *= pow(cf, exp - 1);
268  }
269  }
270 
271  label srRef = 0;
272  rRef = R.rhs()[srRef].index;
273 
274  // Find the matrix element and element position for the rhs
275  pr = kr;
276  for (label s=1; s<Nr; s++)
277  {
278  const label si = R.rhs()[s].index;
279  if (c[si] < c[rRef])
280  {
281  const scalar exp = R.rhs()[srRef].exponent;
282  pr *= pow(max(c[rRef], 0.0), exp);
283  rRef = si;
284  srRef = s;
285  }
286  else
287  {
288  const scalar exp = R.rhs()[s].exponent;
289  pr *= pow(max(c[si], 0.0), exp);
290  }
291  }
292  cr = max(c[rRef], 0.0);
293 
294  {
295  const scalar exp = R.rhs()[srRef].exponent;
296  if (exp < 1)
297  {
298  if (cr > SMALL)
299  {
300  pr *= pow(cr, exp - 1);
301  }
302  else
303  {
304  pr = 0;
305  }
306  }
307  else
308  {
309  pr *= pow(cr, exp - 1);
310  }
311  }
312 
313  return pf*cf - pr*cr;
314 }
315 
316 
317 template<class ReactionThermo, class ThermoType>
319 (
320  const scalar time,
321  const scalarField& c,
322  scalarField& dcdt
323 ) const
324 {
325  const bool reduced = mechRed_->active();
326 
327  const scalar T = c[this->nSpecie_];
328  const scalar p = c[this->nSpecie_ + 1];
329 
330  if (reduced)
331  {
332  // When using DAC, the ODE solver submit a reduced set of species the
333  // complete set is used and only the species in the simplified mechanism
334  // are updated
335  this->c_ = completeC_;
336 
337  // Update the concentration of the species in the simplified mechanism
338  // the other species remain the same and are used only for third-body
339  // efficiencies
340  for (label i=0; i<NsDAC_; i++)
341  {
342  this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
343  }
344  }
345  else
346  {
347  for (label i=0; i<this->nSpecie(); i++)
348  {
349  this->c_[i] = max(c[i], 0.0);
350  }
351  }
352 
353  omega(this->c_, T, p, dcdt);
354 
355  // Constant pressure
356  // dT/dt = ...
357  scalar rho = 0;
358  for (label i=0; i<this->c_.size(); i++)
359  {
360  const scalar W = this->specieThermo_[i].W();
361  rho += W*this->c_[i];
362  }
363 
364  scalar cp = 0;
365  for (label i=0; i<this->c_.size(); i++)
366  {
367  // cp function returns [J/(kmol K)]
368  cp += this->c_[i]*this->specieThermo_[i].cp(p, T);
369  }
370  cp /= rho;
371 
372  // When mechanism reduction is active
373  // dT is computed on the reduced set since dcdt is null
374  // for species not involved in the simplified mechanism
375  scalar dT = 0;
376  for (label i=0; i<this->nSpecie_; i++)
377  {
378  label si;
379  if (reduced)
380  {
381  si = simplifiedToCompleteIndex_[i];
382  }
383  else
384  {
385  si = i;
386  }
387 
388  // ha function returns [J/kmol]
389  const scalar hi = this->specieThermo_[si].ha(p, T);
390  dT += hi*dcdt[i];
391  }
392  dT /= rho*cp;
393 
394  dcdt[this->nSpecie_] = -dT;
395 
396  // dp/dt = ...
397  dcdt[this->nSpecie_ + 1] = 0;
398 }
399 
400 
401 template<class ReactionThermo, class ThermoType>
403 (
404  const scalar t,
405  const scalarField& c,
406  scalarSquareMatrix& dfdc
407 ) const
408 {
409  const bool reduced = mechRed_->active();
410 
411  // If the mechanism reduction is active, the computed Jacobian
412  // is compact (size of the reduced set of species)
413  // but according to the information of the complete set
414  // (i.e. for the third-body efficiencies)
415 
416  const label nSpecie = this->nSpecie_;
417 
418  const scalar T = c[nSpecie];
419  const scalar p = c[nSpecie + 1];
420 
421  if (reduced)
422  {
423  this->c_ = completeC_;
424  for (label i=0; i<NsDAC_; ++i)
425  {
426  this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
427  }
428  }
429  else
430  {
431  forAll(this->c_, i)
432  {
433  this->c_[i] = max(c[i], 0.0);
434  }
435  }
436 
437  dfdc = Zero;
438 
439  forAll(this->reactions_, ri)
440  {
441  if (!reactionsDisabled_[ri])
442  {
443  const Reaction<ThermoType>& R = this->reactions_[ri];
444 
445  const scalar kf0 = R.kf(p, T, this->c_);
446  const scalar kr0 = R.kr(kf0, p, T, this->c_);
447 
448  forAll(R.lhs(), j)
449  {
450  label sj = R.lhs()[j].index;
451  if (reduced)
452  {
453  sj = completeToSimplifiedIndex_[sj];
454  }
455  scalar kf = kf0;
456  forAll(R.lhs(), i)
457  {
458  const label si = R.lhs()[i].index;
459  const scalar el = R.lhs()[i].exponent;
460  if (i == j)
461  {
462  if (el < 1)
463  {
464  if (this->c_[si] > SMALL)
465  {
466  kf *= el*pow(this->c_[si], el - 1);
467  }
468  else
469  {
470  kf = 0;
471  }
472  }
473  else
474  {
475  kf *= el*pow(this->c_[si], el - 1);
476  }
477  }
478  else
479  {
480  kf *= pow(this->c_[si], el);
481  }
482  }
483 
484  forAll(R.lhs(), i)
485  {
486  label si = R.lhs()[i].index;
487  if (reduced)
488  {
489  si = completeToSimplifiedIndex_[si];
490  }
491  const scalar sl = R.lhs()[i].stoichCoeff;
492  dfdc(si, sj) -= sl*kf;
493  }
494  forAll(R.rhs(), i)
495  {
496  label si = R.rhs()[i].index;
497  if (reduced)
498  {
499  si = completeToSimplifiedIndex_[si];
500  }
501  const scalar sr = R.rhs()[i].stoichCoeff;
502  dfdc(si, sj) += sr*kf;
503  }
504  }
505 
506  forAll(R.rhs(), j)
507  {
508  label sj = R.rhs()[j].index;
509  if (reduced)
510  {
511  sj = completeToSimplifiedIndex_[sj];
512  }
513  scalar kr = kr0;
514  forAll(R.rhs(), i)
515  {
516  const label si = R.rhs()[i].index;
517  const scalar er = R.rhs()[i].exponent;
518  if (i == j)
519  {
520  if (er < 1)
521  {
522  if (this->c_[si] > SMALL)
523  {
524  kr *= er*pow(this->c_[si], er - 1);
525  }
526  else
527  {
528  kr = 0;
529  }
530  }
531  else
532  {
533  kr *= er*pow(this->c_[si], er - 1);
534  }
535  }
536  else
537  {
538  kr *= pow(this->c_[si], er);
539  }
540  }
541 
542  forAll(R.lhs(), i)
543  {
544  label si = R.lhs()[i].index;
545  if (reduced)
546  {
547  si = completeToSimplifiedIndex_[si];
548  }
549  const scalar sl = R.lhs()[i].stoichCoeff;
550  dfdc(si, sj) += sl*kr;
551  }
552  forAll(R.rhs(), i)
553  {
554  label si = R.rhs()[i].index;
555  if (reduced)
556  {
557  si = completeToSimplifiedIndex_[si];
558  }
559  const scalar sr = R.rhs()[i].stoichCoeff;
560  dfdc(si, sj) -= sr*kr;
561  }
562  }
563  }
564  }
565 
566  // Calculate the dcdT elements numerically
567  const scalar delta = 1e-3;
568 
569  omega(this->c_, T + delta, p, this->dcdt_);
570  for (label i=0; i<nSpecie; ++i)
571  {
572  dfdc(i, nSpecie) = this->dcdt_[i];
573  }
574 
575  omega(this->c_, T - delta, p, this->dcdt_);
576  for (label i=0; i<nSpecie; ++i)
577  {
578  dfdc(i, nSpecie) = 0.5*(dfdc(i, nSpecie) - this->dcdt_[i])/delta;
579  }
580 
581  dfdc(nSpecie, nSpecie) = 0;
582  dfdc(nSpecie + 1, nSpecie) = 0;
583 }
584 
585 
586 template<class ReactionThermo, class ThermoType>
588 (
589  const scalar t,
590  const scalarField& c,
591  scalarField& dcdt,
592  scalarSquareMatrix& dfdc
593 ) const
594 {
595  jacobian(t, c, dfdc);
596 
597  const scalar T = c[this->nSpecie_];
598  const scalar p = c[this->nSpecie_ + 1];
599 
600  // Note: Uses the c_ field initialized by the call to jacobian above
601  omega(this->c_, T, p, dcdt);
602 }
603 
604 
605 template<class ReactionThermo, class ThermoType>
606 template<class DeltaTType>
608 (
609  const DeltaTType& deltaT
610 )
611 {
612  // Increment counter of time-step
613  timeSteps_++;
614 
615  const bool reduced = mechRed_->active();
616 
617  label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
618 
619  basicSpecieMixture& composition = this->thermo().composition();
620 
621  // CPU time analysis
622  const clockTime clockTime_= clockTime();
623  clockTime_.timeIncrement();
624  scalar reduceMechCpuTime_ = 0;
625  scalar addNewLeafCpuTime_ = 0;
626  scalar growCpuTime_ = 0;
627  scalar solveChemistryCpuTime_ = 0;
628  scalar searchISATCpuTime_ = 0;
629 
630  this->resetTabulationResults();
631 
632  // Average number of active species
633  scalar nActiveSpecies = 0;
634  scalar nAvg = 0;
635 
637 
638  scalar deltaTMin = GREAT;
639 
640  if (!this->chemistry_)
641  {
642  return deltaTMin;
643  }
644 
645  const volScalarField rho
646  (
647  IOobject
648  (
649  "rho",
650  this->time().timeName(),
651  this->mesh(),
652  IOobject::NO_READ,
653  IOobject::NO_WRITE,
654  false
655  ),
656  this->thermo().rho()
657  );
658 
659  const scalarField& T = this->thermo().T();
660  const scalarField& p = this->thermo().p();
661 
662  scalarField c(this->nSpecie_);
663  scalarField c0(this->nSpecie_);
664 
665  // Composition vector (Yi, T, p)
666  scalarField phiq(this->nEqns() + nAdditionalEqn);
667 
668  scalarField Rphiq(this->nEqns() + nAdditionalEqn);
669 
670  forAll(rho, celli)
671  {
672  const scalar rhoi = rho[celli];
673  scalar pi = p[celli];
674  scalar Ti = T[celli];
675 
676  for (label i=0; i<this->nSpecie_; i++)
677  {
678  const volScalarField& Yi = this->Y_[i];
679  c[i] = rhoi*Yi[celli]/this->specieThermo_[i].W();
680  c0[i] = c[i];
681  phiq[i] = Yi[celli];
682  }
683  phiq[this->nSpecie()] = Ti;
684  phiq[this->nSpecie() + 1] = pi;
685  if (tabulation_->variableTimeStep())
686  {
687  phiq[this->nSpecie() + 2] = deltaT[celli];
688  }
689 
690 
691  // Initialise time progress
692  scalar timeLeft = deltaT[celli];
693 
694  // Not sure if this is necessary
695  Rphiq = Zero;
696 
697  clockTime_.timeIncrement();
698 
699  // When tabulation is active (short-circuit evaluation for retrieve)
700  // It first tries to retrieve the solution of the system with the
701  // information stored through the tabulation method
702  if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
703  {
704  // Retrieved solution stored in Rphiq
705  for (label i=0; i<this->nSpecie(); ++i)
706  {
707  c[i] = rhoi*Rphiq[i]/this->specieThermo_[i].W();
708  }
709 
710  searchISATCpuTime_ += clockTime_.timeIncrement();
711  }
712  // This position is reached when tabulation is not used OR
713  // if the solution is not retrieved.
714  // In the latter case, it adds the information to the tabulation
715  // (it will either expand the current data or add a new stored point).
716  else
717  {
718  // Store total time waiting to attribute to add or grow
719  scalar timeTmp = clockTime_.timeIncrement();
720 
721  if (reduced)
722  {
723  // Reduce mechanism change the number of species (only active)
724  mechRed_->reduceMechanism(c, Ti, pi);
725  nActiveSpecies += mechRed_->NsSimp();
726  ++nAvg;
727  scalar timeIncr = clockTime_.timeIncrement();
728  reduceMechCpuTime_ += timeIncr;
729  timeTmp += timeIncr;
730  }
731 
732  // Calculate the chemical source terms
733  while (timeLeft > SMALL)
734  {
735  scalar dt = timeLeft;
736  if (reduced)
737  {
738  // completeC_ used in the overridden ODE methods
739  // to update only the active species
740  completeC_ = c;
741 
742  // Solve the reduced set of ODE
743  this->solve
744  (
745  simplifiedC_, Ti, pi, dt, this->deltaTChem_[celli]
746  );
747 
748  for (label i=0; i<NsDAC_; ++i)
749  {
750  c[simplifiedToCompleteIndex_[i]] = simplifiedC_[i];
751  }
752  }
753  else
754  {
755  this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
756  }
757  timeLeft -= dt;
758  }
759 
760  {
761  scalar timeIncr = clockTime_.timeIncrement();
762  solveChemistryCpuTime_ += timeIncr;
763  timeTmp += timeIncr;
764  }
765 
766  // If tabulation is used, we add the information computed here to
767  // the stored points (either expand or add)
768  if (tabulation_->active())
769  {
770  forAll(c, i)
771  {
772  Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W();
773  }
774  if (tabulation_->variableTimeStep())
775  {
776  Rphiq[Rphiq.size()-3] = Ti;
777  Rphiq[Rphiq.size()-2] = pi;
778  Rphiq[Rphiq.size()-1] = deltaT[celli];
779  }
780  else
781  {
782  Rphiq[Rphiq.size()-2] = Ti;
783  Rphiq[Rphiq.size()-1] = pi;
784  }
785  label growOrAdd =
786  tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]);
787 
788  if (growOrAdd)
789  {
790  this->setTabulationResultsAdd(celli);
791  addNewLeafCpuTime_ += clockTime_.timeIncrement() + timeTmp;
792  }
793  else
794  {
795  this->setTabulationResultsGrow(celli);
796  growCpuTime_ += clockTime_.timeIncrement() + timeTmp;
797  }
798  }
799 
800  // When operations are done and if mechanism reduction is active,
801  // the number of species (which also affects nEqns) is set back
802  // to the total number of species (stored in the mechRed object)
803  if (reduced)
804  {
805  this->nSpecie_ = mechRed_->nSpecie();
806  }
807  deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
808 
809  this->deltaTChem_[celli] =
810  min(this->deltaTChem_[celli], this->deltaTChemMax_);
811  }
812 
813  // Set the RR vector (used in the solver)
814  for (label i=0; i<this->nSpecie_; ++i)
815  {
816  this->RR_[i][celli] =
817  (c[i] - c0[i])*this->specieThermo_[i].W()/deltaT[celli];
818  }
819  }
820 
821  if (mechRed_->log() || tabulation_->log())
822  {
823  cpuSolveFile_()
824  << this->time().timeOutputValue()
825  << " " << solveChemistryCpuTime_ << endl;
826  }
827 
828  if (mechRed_->log())
829  {
830  cpuReduceFile_()
831  << this->time().timeOutputValue()
832  << " " << reduceMechCpuTime_ << endl;
833  }
834 
835  if (tabulation_->active())
836  {
837  // Every time-step, look if the tabulation should be updated
838  tabulation_->update();
839 
840  // Write the performance of the tabulation
841  tabulation_->writePerformance();
842 
843  if (tabulation_->log())
844  {
845  cpuRetrieveFile_()
846  << this->time().timeOutputValue()
847  << " " << searchISATCpuTime_ << endl;
848 
849  cpuGrowFile_()
850  << this->time().timeOutputValue()
851  << " " << growCpuTime_ << endl;
852 
853  cpuAddFile_()
854  << this->time().timeOutputValue()
855  << " " << addNewLeafCpuTime_ << endl;
856  }
857  }
858 
859  if (reduced && nAvg && mechRed_->log())
860  {
861  // Write average number of species
862  nActiveSpeciesFile_()
863  << this->time().timeOutputValue()
864  << " " << nActiveSpecies/nAvg << endl;
865  }
866 
867  if (Pstream::parRun())
868  {
869  List<bool> active(composition.active());
870  Pstream::listCombineGather(active, orEqOp<bool>());
871  Pstream::listCombineScatter(active);
872 
873  forAll(active, i)
874  {
875  if (active[i])
876  {
877  composition.setActive(i);
878  }
879  }
880  }
881 
882  forAll(this->Y(), i)
883  {
884  if (composition.active(i))
885  {
886  this->Y()[i].writeOpt(IOobject::AUTO_WRITE);
887  }
888  }
889 
890  return deltaTMin;
891 }
892 
893 
894 template<class ReactionThermo, class ThermoType>
896 (
897  const scalar deltaT
898 )
899 {
900  // Don't allow the time-step to change more than a factor of 2
901  return min
902  (
904  2*deltaT
905  );
906 }
907 
908 
909 template<class ReactionThermo, class ThermoType>
911 (
912  const scalarField& deltaT
913 )
914 {
915  return this->solve<scalarField>(deltaT);
916 }
917 
918 
919 template<class ReactionThermo, class ThermoType>
922 (
923  const label celli
924 )
925 {
926  tabulationResults_[celli] = 0.0;
927 }
928 
929 
930 template<class ReactionThermo, class ThermoType>
932 setTabulationResultsGrow(const label celli)
933 {
934  tabulationResults_[celli] = 1.0;
935 }
936 
937 
938 template<class ReactionThermo, class ThermoType>
941 (
942  const label celli
943 )
944 {
945  tabulationResults_[celli] = 2.0;
946 }
947 
948 
949 // ************************************************************************* //
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::TDACChemistryModel::~TDACChemistryModel
virtual ~TDACChemistryModel()
Destructor.
Definition: TDACChemistryModel.C:148
Foam::TDACChemistryModel::derivatives
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Definition: TDACChemistryModel.C:319
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))
Definition: gmvOutputSpray.H:25
cp
const volScalarField & cp
Definition: setRegionSolidFields.H:8
Foam::basicSpecieMixture
Specialization of basicMultiComponentMixture for a mixture consisting of a number for molecular speci...
Definition: basicSpecieMixture.H:58
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::IOobject::typeHeaderOk
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
Definition: IOobjectTemplates.C:39
Foam::clockTime::timeIncrement
double timeIncrement() const
The time [seconds] since the last call to timeIncrement()
Definition: clockTimeI.H:60
TDACChemistryModel.H
thermo
psiReactionThermo & thermo
Definition: createFields.H:28
Foam::reactingMixture
Foam::reactingMixture.
Definition: reactingMixture.H:57
nSpecie
label nSpecie
Definition: readInitialConditions.H:17
thermo
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
Foam::orEqOp
Definition: ops.H:86
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::TDACChemistryModel::active
bool active(const label i) const
Definition: TDACChemistryModelI.H:91
Foam::exp
dimensionedScalar exp(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:261
rho
rho
Definition: readInitialConditions.H:88
composition
basicSpecieMixture & composition
Definition: createFieldRefs.H:6
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
solve
CEqn solve()
Foam::TDACChemistryModel::setTabulationResultsGrow
void setTabulationResultsGrow(const label celli)
Definition: TDACChemistryModel.C:932
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::UniformField
A class representing the concept of a uniform field which stores only the single value and providing ...
Definition: UniformField.H:51
R
#define R(A, B, C, D, E, F, K, M)
localEulerDdtScheme.H
Foam::Field< scalar >
Foam::hashedWordList
A wordList with hashed named lookup, which can be faster in some situations than using the normal lis...
Definition: hashedWordList.H:54
controlDict
runTime controlDict().readEntry("adjustTimeStep"
Definition: debug.C:143
correct
fvOptions correct(rho)
delta
scalar delta
Definition: LISASMDCalcMethod2.H:8
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:42
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
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
timeName
word timeName
Definition: getTimeIndex.H:3
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::TDACChemistryModel::omega
virtual void omega(const scalarField &c, const scalar T, const scalar p, scalarField &dcdt) const
dc/dt = omega, rate of change in concentration, for each species
Definition: TDACChemistryModel.C:156
T
const volScalarField & T
Definition: createFieldRefs.H:2
Y
PtrList< volScalarField > & Y
Definition: createFieldRefs.H:7
Foam::New
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Definition: DimensionedFieldReuseFunctions.H:105
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::SquareMatrix< scalar >
Foam::TDACChemistryModel::setTabulationResultsAdd
void setTabulationResultsAdd(const label celli)
Definition: TDACChemistryModel.C:922
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::autoPtr::ref
T & ref()
Return reference to the managed object without nullptr checking.
Definition: autoPtr.H:161
Foam::dynamicCast
To & dynamicCast(From &r)
Definition: typeInfo.H:88
Foam::List< bool >
Foam::TDACChemistryModel::jacobian
void jacobian(const scalar t, const scalarField &c, scalarSquareMatrix &dfdc) const
Pure jacobian function for tabulation.
Definition: TDACChemistryModel.C:403
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::clockTime
Starts timing and returns elapsed time from start. Uses std::chrono::high_resolution_clock for better...
Definition: clockTime.H:62
Foam::TDACChemistryModel
Extends StandardChemistryModel by adding the TDAC method.
Definition: chemistryReductionMethod.H:49
Foam::GeometricField< scalar, fvPatchField, volMesh >
Foam::Reaction
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:59
clockTime.H
UniformField.H
Foam::dimless
const dimensionSet dimless
Dimensionless.
Definition: dimensionSets.C:189
Foam::StandardChemistryModel
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
Definition: StandardChemistryModel.H:61
Foam::TDACChemistryModel::setTabulationResultsRetrieve
void setTabulationResultsRetrieve(const label celli)
Definition: TDACChemistryModel.C:941