Skip to content

Commit 1012ad1

Browse files
OpenFOAMFansOpenFOAMFans
authored
Merge pull request #268 from wengzf20/canteraErr
catch error from Cantera and report in OpenFOAM
2 parents 217c472 + e877355 commit 1012ad1

1 file changed

Lines changed: 127 additions & 118 deletions

File tree

‎src/dfChemistryModel/dfChemistryModel.C‎

Lines changed: 127 additions & 118 deletions
Original file line numberDiff line numberDiff line change
@@ -373,181 +373,190 @@ void Foam::dfChemistryModel<ThermoType>::setNumerics(Cantera::ReactorNet &sim)
373373
template<class ThermoType>
374374
void Foam::dfChemistryModel<ThermoType>::correctThermo()
375375
{
376-
psi_.oldTime();
377-
378-
forAll(T_, celli)
376+
try
379377
{
380-
forAll(Y_, i)
381-
{
382-
yTemp_[i] = Y_[i][celli];
383-
}
384-
CanteraGas_->setState_PY(p_[celli], yTemp_.begin());
385-
if(mixture_.heName()=="ha")
386-
{
387-
CanteraGas_->setState_HP(thermo_.he()[celli], p_[celli]); // setState_HP needs (J/kg)
388-
}
389-
else if(mixture_.heName()=="ea")
378+
psi_.oldTime();
379+
380+
forAll(T_, celli)
390381
{
391-
scalar ha = thermo_.he()[celli] + p_[celli]/rho_[celli];
392-
CanteraGas_->setState_HP(ha, p_[celli]);
393-
}
382+
forAll(Y_, i)
383+
{
384+
yTemp_[i] = Y_[i][celli];
385+
}
386+
CanteraGas_->setState_PY(p_[celli], yTemp_.begin());
387+
if(mixture_.heName()=="ha")
388+
{
389+
CanteraGas_->setState_HP(thermo_.he()[celli], p_[celli]); // setState_HP needs (J/kg)
390+
}
391+
else if(mixture_.heName()=="ea")
392+
{
393+
scalar ha = thermo_.he()[celli] + p_[celli]/rho_[celli];
394+
CanteraGas_->setState_HP(ha, p_[celli]);
395+
}
394396

395397

396-
T_[celli] = CanteraGas_->temperature();
398+
T_[celli] = CanteraGas_->temperature();
397399

398-
psi_[celli] = mixture_.psi(p_[celli],T_[celli]);
400+
psi_[celli] = mixture_.psi(p_[celli],T_[celli]);
399401

400-
rho_[celli] = mixture_.rho(p_[celli],T_[celli]);
402+
rho_[celli] = mixture_.rho(p_[celli],T_[celli]);
401403

402-
mu_[celli] = mixture_.CanteraTransport()->viscosity(); // Pa-s
404+
mu_[celli] = mixture_.CanteraTransport()->viscosity(); // Pa-s
403405

404-
alpha_[celli] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass()); // kg/(m*s)
405-
// thermalConductivity() W/m/K
406-
// cp_mass() J/kg/K
406+
alpha_[celli] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass()); // kg/(m*s)
407+
// thermalConductivity() W/m/K
408+
// cp_mass() J/kg/K
407409

408-
if (mixture_.transportModelName() == "UnityLewis")
409-
{
410-
forAll(rhoD_, i)
410+
if (mixture_.transportModelName() == "UnityLewis")
411411
{
412-
rhoD_[i][celli] = alpha_[celli];
412+
forAll(rhoD_, i)
413+
{
414+
rhoD_[i][celli] = alpha_[celli];
415+
}
413416
}
414-
}
415-
else
416-
{
417-
mixture_.CanteraTransport()->getMixDiffCoeffsMass(dTemp_.begin()); // m2/s
418-
419-
CanteraGas_->getEnthalpy_RT(hrtTemp_.begin()); //hrtTemp_=m_h0_RT non-dimension
420-
// constant::physicoChemical::R.value() J/(mol·k)
421-
const scalar RT = constant::physicoChemical::R.value()*1e3*T_[celli]; // J/kmol/K
422-
forAll(rhoD_, i)
417+
else
423418
{
424-
rhoD_[i][celli] = rho_[celli]*dTemp_[i];
419+
mixture_.CanteraTransport()->getMixDiffCoeffsMass(dTemp_.begin()); // m2/s
425420

426-
// CanteraGas_->molecularWeight(i) kg/kmol
427-
hai_[i][celli] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i);
421+
CanteraGas_->getEnthalpy_RT(hrtTemp_.begin()); //hrtTemp_=m_h0_RT non-dimension
422+
// constant::physicoChemical::R.value() J/(mol·k)
423+
const scalar RT = constant::physicoChemical::R.value()*1e3*T_[celli]; // J/kmol/K
424+
forAll(rhoD_, i)
425+
{
426+
rhoD_[i][celli] = rho_[celli]*dTemp_[i];
427+
428+
// CanteraGas_->molecularWeight(i) kg/kmol
429+
hai_[i][celli] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i);
430+
}
428431
}
429432
}
430-
}
431433

432434

433-
const volScalarField::Boundary& pBf = p_.boundaryField();
435+
const volScalarField::Boundary& pBf = p_.boundaryField();
434436

435-
volScalarField::Boundary& rhoBf = rho_.boundaryFieldRef();
437+
volScalarField::Boundary& rhoBf = rho_.boundaryFieldRef();
436438

437-
volScalarField::Boundary& TBf = T_.boundaryFieldRef();
439+
volScalarField::Boundary& TBf = T_.boundaryFieldRef();
438440

439-
volScalarField::Boundary& psiBf = psi_.boundaryFieldRef();
441+
volScalarField::Boundary& psiBf = psi_.boundaryFieldRef();
440442

441-
volScalarField::Boundary& hBf = thermo_.he().boundaryFieldRef();
443+
volScalarField::Boundary& hBf = thermo_.he().boundaryFieldRef();
442444

443-
volScalarField::Boundary& muBf = mu_.boundaryFieldRef();
445+
volScalarField::Boundary& muBf = mu_.boundaryFieldRef();
444446

445-
volScalarField::Boundary& alphaBf = alpha_.boundaryFieldRef();
447+
volScalarField::Boundary& alphaBf = alpha_.boundaryFieldRef();
446448

447-
forAll(T_.boundaryField(), patchi)
448-
{
449-
const fvPatchScalarField& pp = pBf[patchi];
450-
fvPatchScalarField& prho = rhoBf[patchi];
451-
fvPatchScalarField& pT = TBf[patchi];
452-
fvPatchScalarField& ppsi = psiBf[patchi];
453-
fvPatchScalarField& ph = hBf[patchi];
454-
fvPatchScalarField& pmu = muBf[patchi];
455-
fvPatchScalarField& palpha = alphaBf[patchi];
456-
457-
if (pT.fixesValue())
449+
forAll(T_.boundaryField(), patchi)
458450
{
459-
forAll(pT, facei)
451+
const fvPatchScalarField& pp = pBf[patchi];
452+
fvPatchScalarField& prho = rhoBf[patchi];
453+
fvPatchScalarField& pT = TBf[patchi];
454+
fvPatchScalarField& ppsi = psiBf[patchi];
455+
fvPatchScalarField& ph = hBf[patchi];
456+
fvPatchScalarField& pmu = muBf[patchi];
457+
fvPatchScalarField& palpha = alphaBf[patchi];
458+
459+
if (pT.fixesValue())
460460
{
461-
forAll(Y_, i)
461+
forAll(pT, facei)
462462
{
463-
yTemp_[i] = Y_[i].boundaryField()[patchi][facei];
464-
}
465-
CanteraGas_->setState_TPY(pT[facei], pp[facei], yTemp_.begin());
463+
forAll(Y_, i)
464+
{
465+
yTemp_[i] = Y_[i].boundaryField()[patchi][facei];
466+
}
467+
CanteraGas_->setState_TPY(pT[facei], pp[facei], yTemp_.begin());
466468

467-
ph[facei] = CanteraGas_->enthalpy_mass();
469+
ph[facei] = CanteraGas_->enthalpy_mass();
468470

469-
ppsi[facei] = mixture_.psi(pp[facei],pT[facei]);
471+
ppsi[facei] = mixture_.psi(pp[facei],pT[facei]);
470472

471-
prho[facei] = mixture_.rho(pp[facei],pT[facei]);
473+
prho[facei] = mixture_.rho(pp[facei],pT[facei]);
472474

473-
pmu[facei] = mixture_.CanteraTransport()->viscosity();
475+
pmu[facei] = mixture_.CanteraTransport()->viscosity();
474476

475-
palpha[facei] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass());
477+
palpha[facei] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass());
476478

477-
if (mixture_.transportModelName() == "UnityLewis")
478-
{
479-
forAll(rhoD_, i)
479+
if (mixture_.transportModelName() == "UnityLewis")
480480
{
481-
rhoD_[i].boundaryFieldRef()[patchi][facei] = palpha[facei];
481+
forAll(rhoD_, i)
482+
{
483+
rhoD_[i].boundaryFieldRef()[patchi][facei] = palpha[facei];
484+
}
482485
}
483-
}
484-
else
485-
{
486-
mixture_.CanteraTransport()->getMixDiffCoeffsMass(dTemp_.begin());
487-
488-
CanteraGas_->getEnthalpy_RT(hrtTemp_.begin());
489-
const scalar RT = constant::physicoChemical::R.value()*1e3*pT[facei];
490-
forAll(rhoD_, i)
486+
else
491487
{
492-
rhoD_[i].boundaryFieldRef()[patchi][facei] = prho[facei]*dTemp_[i];
488+
mixture_.CanteraTransport()->getMixDiffCoeffsMass(dTemp_.begin());
489+
490+
CanteraGas_->getEnthalpy_RT(hrtTemp_.begin());
491+
const scalar RT = constant::physicoChemical::R.value()*1e3*pT[facei];
492+
forAll(rhoD_, i)
493+
{
494+
rhoD_[i].boundaryFieldRef()[patchi][facei] = prho[facei]*dTemp_[i];
493495

494-
hai_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i);
496+
hai_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i);
497+
}
495498
}
496499
}
497500
}
498-
}
499-
else
500-
{
501-
forAll(pT, facei)
501+
else
502502
{
503-
forAll(Y_, i)
504-
{
505-
yTemp_[i] = Y_[i].boundaryField()[patchi][facei];
506-
}
507-
CanteraGas_->setState_PY(pp[facei], yTemp_.begin());
508-
if(mixture_.heName()=="ha")
509-
{
510-
CanteraGas_->setState_HP(ph[facei], pp[facei]);
511-
}
512-
else if(mixture_.heName()=="ea")
503+
forAll(pT, facei)
513504
{
514-
scalar ha = ph[facei] + pp[facei]/prho[facei];
515-
CanteraGas_->setState_HP(ha, pp[facei]);
516-
}
505+
forAll(Y_, i)
506+
{
507+
yTemp_[i] = Y_[i].boundaryField()[patchi][facei];
508+
}
509+
CanteraGas_->setState_PY(pp[facei], yTemp_.begin());
510+
if(mixture_.heName()=="ha")
511+
{
512+
CanteraGas_->setState_HP(ph[facei], pp[facei]);
513+
}
514+
else if(mixture_.heName()=="ea")
515+
{
516+
scalar ha = ph[facei] + pp[facei]/prho[facei];
517+
CanteraGas_->setState_HP(ha, pp[facei]);
518+
}
517519

518-
pT[facei] = CanteraGas_->temperature();
520+
pT[facei] = CanteraGas_->temperature();
519521

520-
ppsi[facei] = mixture_.psi(pp[facei],pT[facei]);
522+
ppsi[facei] = mixture_.psi(pp[facei],pT[facei]);
521523

522-
prho[facei] = mixture_.rho(pp[facei],pT[facei]);
524+
prho[facei] = mixture_.rho(pp[facei],pT[facei]);
523525

524-
pmu[facei] = mixture_.CanteraTransport()->viscosity();
526+
pmu[facei] = mixture_.CanteraTransport()->viscosity();
525527

526-
palpha[facei] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass());
528+
palpha[facei] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass());
527529

528-
if (mixture_.transportModelName() == "UnityLewis")
529-
{
530-
forAll(rhoD_, i)
530+
if (mixture_.transportModelName() == "UnityLewis")
531531
{
532-
rhoD_[i].boundaryFieldRef()[patchi][facei] = palpha[facei];
532+
forAll(rhoD_, i)
533+
{
534+
rhoD_[i].boundaryFieldRef()[patchi][facei] = palpha[facei];
535+
}
533536
}
534-
}
535-
else
536-
{
537-
mixture_.CanteraTransport()->getMixDiffCoeffsMass(dTemp_.begin());
538-
539-
CanteraGas_->getEnthalpy_RT(hrtTemp_.begin());
540-
const scalar RT = constant::physicoChemical::R.value()*1e3*pT[facei];
541-
forAll(rhoD_, i)
537+
else
542538
{
543-
rhoD_[i].boundaryFieldRef()[patchi][facei] = prho[facei]*dTemp_[i];
539+
mixture_.CanteraTransport()->getMixDiffCoeffsMass(dTemp_.begin());
544540

545-
hai_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i);
541+
CanteraGas_->getEnthalpy_RT(hrtTemp_.begin());
542+
const scalar RT = constant::physicoChemical::R.value()*1e3*pT[facei];
543+
forAll(rhoD_, i)
544+
{
545+
rhoD_[i].boundaryFieldRef()[patchi][facei] = prho[facei]*dTemp_[i];
546+
547+
hai_[i].boundaryFieldRef()[patchi][facei] = hrtTemp_[i]*RT/CanteraGas_->molecularWeight(i);
548+
}
546549
}
547550
}
548551
}
549552
}
550553
}
554+
catch(Cantera::CanteraError& err)
555+
{
556+
std::cerr << err.what() << '\n';
557+
FatalErrorInFunction
558+
<< abort(FatalError);
559+
}
551560
}
552561

553562
template<class ThermoType>

0 commit comments

Comments
 (0)