print thermo data and reaction data to stdout

This commit is contained in:
Yeongdo Park 2018-11-08 13:58:29 -05:00
parent 8e31407c61
commit 3185825983
2 changed files with 249 additions and 215 deletions

377
ct2foam.C
View file

@ -34,6 +34,8 @@ Description
#include <cantera/transport.h>
#include <cantera/IdealGasMix.h>
#include "ofFormats.h"
#include "fvCFD.H"
#include "psiReactionThermo.H"
#include "psiChemistryModel.H"
@ -45,16 +47,16 @@ Description
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
const scalar relTol = 10.0 / 100.;
inline scalar relError (scalar x, scalar x0)
std::string stringNASACoefs (doublereal *coef)
{
return (x - x0) / x0;
}
std::string arr;
inline bool exceedTolerence (scalar rError)
{
return rError > relTol || rError < -relTol;
for (label j = 0; j < 7; j++)
{
arr += fmt::format(" {:15.10E} ", coef[j]);
}
return arr;
}
int main(int argc, char *argv[])
@ -66,265 +68,210 @@ int main(int argc, char *argv[])
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCreate Cantera object\n" << endl;
Cantera::IdealGasMix gas_ ("gri30.cti", "gri30_mix");
Cantera::Transport *tr_ = Cantera::newTransportMgr("Mix", &gas_);
label nCanteraSp_ = gas_.nSpecies();
int nsp = gas_.nSpecies();
int type;
doublereal c[15];
doublereal minTemp, maxTemp, refPressure;
Cantera::MultiSpeciesThermo& sp = gas_.speciesThermo();
int n, j;
// these constants define the location of coefficient "a6" in the
// cofficient array c. The c array contains Tmid in the first
// location, followed by the 7 low-temperature coefficients, then
// the seven high-temperature ones.
for (n = 0; n < nsp; n++)
for (size_t n = 0; n < gas_.nSpecies(); n++)
{
Cantera::writelog("\n\n {} (original):", gas_.speciesName(n));
int type;
doublereal c[15];
doublereal minTemp, maxTemp, refPressure;
// get the NASA coefficients in array c
sp.reportParams(n, type, c, minTemp, maxTemp, refPressure);
gas_.speciesThermo().reportParams(n, type, c, minTemp, maxTemp, refPressure);
// print the unmodified NASA coefficients
Cantera::writelog("\n ");
for (j = 1; j < 8; j++) {
Cantera::writelog(" A{} ", j);
}
Cantera::writelog("\n low:");
for (j = 1; j < 8; j++) {
Cantera::writelog(" {:10.4E} ", c[j]);
}
Cantera::writelog("\n high:");
for (j = 8; j < 15; j++) {
Cantera::writelog(" {:10.4E} ", c[j]);
}
Cantera::writelog("\n ");
// print the modified NASA coefficients
Cantera::writelog("\n\n {} (modified):", gas_.speciesName(n).c_str());
Cantera::writelog("\n ");
for (j = 1; j < 8; j++) {
Cantera::writelog(" A{} ", j);
}
Cantera::writelog("\n low:");
for (j = 1; j < 8; j++) {
Cantera::writelog(" {:10.4E} ", c[j]);
}
Cantera::writelog("\n high:");
for (j = 8; j < 15; j++) {
Cantera::writelog(" {:10.4E} ", c[j]);
}
Cantera::writelog("\n");
Cantera::writelog(
fmt::format(
thermoFormat,
gas_.speciesName(n),
gas_.molecularWeight(n),
gas_.charge(n),
minTemp,
maxTemp,
c[0],
stringNASACoefs(&c[1]),
stringNASACoefs(&c[8])
)
);
}
label nrxn = gas_.nReactions();
// Write species name list
Cantera::writelog("species\n");
Cantera::writelog("{}\n", gas_.nSpecies());
Cantera::writelog("(\n");
for (size_t n = 0; n < gas_.nSpecies(); n++)
{
Cantera::writelog("{}\n", gas_.speciesName(n));
}
Cantera::writelog(")\n");
Cantera::writelog(";\n\n");
for (label k = 0; k < nrxn; k++)
// Write reaction list
Cantera::writelog( "reactions\n{\n");
for (size_t k = 0; k < gas_.nReactions(); k++)
{
std::shared_ptr<Cantera::Reaction> r(gas_.reaction(k));
std::string irn( r->reversible ? "reversible" : "irreversible" );
std::string rate("Arrhenius");
std::string rxn("Reaction");
std::string ffn;
double c[5] = {0};
// reaction name
Cantera::writelog( fmt::format( " un-named-reaction-{}\n", k));
Cantera::writelog( " {\n" );
switch (gas_.reactionType(k))
{
case Cantera::ELEMENTARY_RXN:
Info << "ELEMENTARY_RXN" << endl;
Info << std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.preExponentialFactor() << tab
<< std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.temperatureExponent() << tab
<< std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.activationEnergy_R() << endl;
Cantera::writelog( fmt::format( rxnTypeFormat, irn + rate + rxn));
Cantera::writelog( fmt::format( rxnEqnFormat,
r->Reaction::reactantString() + " = " + r->Reaction::productString()));
Cantera::writelog(
fmt::format(
arrheniusFormat,
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.preExponentialFactor(),
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.temperatureExponent(),
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.activationEnergy_R()
)
);
break;
case Cantera::THREE_BODY_RXN:
Info << "THREE_BODY_RXN" << endl;
Info << std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.preExponentialFactor() << tab
<< std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.temperatureExponent() << tab
<< std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.activationEnergy_R() << endl;
rate = "thirdBody" + rate;
for (label l = 0; l < nsp; l++)
Cantera::writelog( fmt::format( rxnTypeFormat, irn + rate + rxn));
Cantera::writelog( fmt::format( rxnEqnFormat,
r->Reaction::reactantString() + " = " + r->Reaction::productString()));
Cantera::writelog(
fmt::format(
arrheniusFormat,
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.preExponentialFactor(),
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.temperatureExponent(),
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.activationEnergy_R()
)
);
Cantera::writelog(" coeffs \n");
Cantera::writelog(fmt::format("{}\n(\n", gas_.nSpecies()));
for (size_t l = 0; l < gas_.nSpecies(); l++)
{
Info << gas_.speciesName(l) << tab
<< std::dynamic_pointer_cast<Cantera::ThreeBodyReaction>(r)->third_body.efficiency(gas_.speciesName(l)) << endl;
Cantera::writelog(
fmt::format(
"({} {})\n",
gas_.speciesName(l),
std::dynamic_pointer_cast<Cantera::ThreeBodyReaction>(r)->third_body.efficiency(gas_.speciesName(l))
)
);
}
Cantera::writelog(")\n;\n");
break;
case Cantera::FALLOFF_RXN:
Info << "FALLOFF_RXN" << endl;
case Cantera::CHEMACT_RXN:
Info << std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.preExponentialFactor() << tab;
Info << std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.temperatureExponent() << tab;
Info << std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.activationEnergy_R() << endl;
rxn = (r->reaction_type == Cantera::FALLOFF_RXN ? "FallOff" : "ChemicallyActivated") + rxn;
Info << std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.preExponentialFactor() << tab;
Info << std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.temperatureExponent() << tab;
Info << std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.activationEnergy_R() << endl;
switch (std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getType())
{
case Cantera::SIMPLE_FALLOFF:
rxn = "Lindemann" + rxn;
ffn = lindemannFormat;
break;
Info << std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getType() << tab
<< std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->nParameters() << endl;
case Cantera::TROE_FALLOFF:
rxn = "Troe" + rxn;
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getParameters(c);
ffn = fmt::format(troeFormat, c[0], c[1], c[2], c[3]);
break;
case Cantera::SRI_FALLOFF:
rxn = "SRI" + rxn;
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getParameters(c);
ffn = fmt::format(sriFormat, c[0], c[1], c[2], c[3], c[4]);
break;
}
Cantera::writelog( fmt::format( rxnTypeFormat, irn + rate + rxn));
Cantera::writelog( fmt::format( rxnEqnFormat,
r->Reaction::reactantString() + " = " + r->Reaction::productString()));
Cantera::writelog(
fmt::format(
falloffFormat,
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.preExponentialFactor(),
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.temperatureExponent(),
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.activationEnergy_R(),
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.preExponentialFactor(),
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.temperatureExponent(),
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.activationEnergy_R()
)
);
Cantera::writelog( ffn );
Cantera::writelog( std::string() +
" thirdBodyEfficiencies\n" +
" {\n" +
" coeffs \n"
);
Cantera::writelog(fmt::format("{}\n(\n", gas_.nSpecies()));
for (size_t l = 0; l < gas_.nSpecies(); l++)
{
Cantera::writelog(
fmt::format(
"({} {})\n",
gas_.speciesName(l),
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->third_body.efficiency(gas_.speciesName(l))
)
);
}
Cantera::writelog(")\n;\n");
Cantera::writelog(" }\n");
break;
case Cantera::PLOG_RXN:
Info << "PLOG_RXN" << endl;
break;
case Cantera::CHEBYSHEV_RXN:
Info << "CHEBYSHEV_RXN" << endl;
break;
case Cantera::CHEMACT_RXN:
Info << "CHEMACT_RXN" << endl;
break;
default:
break;
}
Info << gas_.reactionString(k) << endl;
Cantera::writelog( " }\n" );
}
Cantera::writelog( "}\n");
/*
Info<< "\nCreate temperature space\n" << endl;
scalar Tl = 300;
scalar Tu = 3000;
scalarField T(100, 0.0);
forAll (T, i)
{
T[i] = i * (Tu - Tl) / scalar(T.size()) + Tl;
}
// initialize composition array
scalarField XY(nCanteraSp_, 0.0);
forAll(thermo.composition().species(), i)
{
XY[i] = Y[i][0];
}
Info<< "\nTemperature Loop\n" << endl;
label errorCount = 0;
label testCount = 0;
forAll(T, i)
{
// OpenFOAM diffusivity model
scalar h0 = 0.0;
forAll(Y, k)
{
h0 += XY[k]*specieData[k].Hs(p[0], T[i]);
}
thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0);
thermo.correct();
diff.correct();
// Cantera gas transport
gas_.setState_TPY(T[i], p[0], XY.data());
scalarField Dc (composition.species().size(), 0.0);
scalarField Dij (composition.species().size()*composition.species().size(), 0.0);
const scalar relTol = 1.0e-15;
// check for singular case (Xi = 1)
label idxFracOne = -1;
scalar sumX = sum(XY);
forAll (XY, k)
{
if ((sumX-XY[k])/sumX < relTol)
{
idxFracOne = k;
break;
}
}
if (idxFracOne < 0)
{
tr_->getMixDiffCoeffsMass(Dc.data());
}
else
{
tr_->getBinaryDiffCoeffs(nCanteraSp_, Dij.data());
forAll(Dc, k)
{
Dc[k] = Dij[nCanteraSp_*idxFracOne+k];
}
}
forAll(thermo.composition().species(), k)
{
if (exceedTolerence(relError(diff.D(k)[0], Dc[k])))
{
Info << thermo.composition().species()[k] << ", T = " << T[i]
<< " relative error = " << (100. * relError(diff.D(k)[0], Dc[k])) << " %" << endl;
errorCount++;
}
testCount++;
post<< thermo.T()[0] << token::TAB
<< thermo.composition().species()[k] << token::TAB
<< diff.D(k)[0] << token::TAB
<< Dc[k] << token::TAB
<< 100*(relError(diff.D(k)[0], Dc[k])) << endl;
}
if (exceedTolerence(relError(diff.mu()[0], tr_->viscosity())))
{
Info << "mu, T = " << T[i] << " relative error = " << (100. * relError(diff.mu()[0], tr_->viscosity())) << " %"<< endl;
errorCount++;
}
testCount++;
post<< thermo.T()[0] << token::TAB
<< "mu" << token::TAB
<< diff.mu()[0] << token::TAB
<< tr_->viscosity() << token::TAB
<< 100*(relError(diff.mu()[0], tr_->viscosity())) << endl;
if (exceedTolerence(relError(diff.k()[0], tr_->thermalConductivity())))
{
Info << "k, T = " << T[i] << " relative error = " << (100. * relError(diff.k()[0], tr_->thermalConductivity())) << " %"<< endl;
errorCount++;
}
testCount++;
post<< thermo.T()[0] << token::TAB
<< "k" << token::TAB
<< diff.k()[0] << token::TAB
<< tr_->thermalConductivity() << token::TAB
<< 100*(relError(diff.k()[0], tr_->thermalConductivity())) << endl;
}
Info << errorCount << " / " << testCount << " End" << nl << endl;
if (errorCount == 0)
{
return 0;
}
else
{
return -1;
}
*/
return 0;
}

87
ofFormats.h Normal file
View file

@ -0,0 +1,87 @@
#ifndef OF_FORMATS_H
#define OF_FORMATS_H
#include <string>
std::string thermoFormat = std::string() +
"{}\n" + // 0 name
"{{\n" +
" specie\n" +
" {{\n" +
" nMoles 1;\n" +
" molWeight {};\n" + // 1 molecular weight
" nCharges {};\n" + // 2 number of elementary charges
" }}\n" +
" thermodynamics\n" +
" {{\n" +
" Tlow {};\n" + // 3 Cp fit low temperature end
" Thigh {};\n" + // 4 Cp fit high temperature end
" Tcommon {};\n" + // 5 Cp fit intermediate temperature
" highCpCoeffs ( {} );\n" + // 6 Cp fit coefficients in high temperature range
" lowCpCoeffs ( {} );\n" + // 7 Cp fit coefficients in low temperature range
" }}\n" +
" transport\n" +
" {{\n" +
" As 1.67212e-06;\n" +
" Ts 170.672;\n" +
" }}\n" +
"}}\n"
;
std::string arrheniusFormat = std::string() +
" A {};\n" +
" beta {};\n" +
" Ta {};\n"
;
std::string falloffFormat = std::string() +
" k0\n" +
" {{\n" +
" A {};\n" +
" beta {};\n" +
" Ta {};\n" +
" }}\n" +
" kInf\n" +
" {{\n" +
" A {};\n" +
" beta {};\n" +
" Ta {};\n" +
" }}\n"
;
std::string lindemannFormat = std::string() +
" F\n" +
" {\n" +
" }\n"
;
std::string troeFormat = std::string() +
" F\n" +
" {{\n" +
" alpha {};\n" +
" Tsss {};\n" +
" Ts {};\n" +
" Tss {};\n" +
" }}\n"
;
std::string sriFormat = std::string() +
" F\n" +
" {{\n" +
" a {};\n" +
" b {};\n" +
" c {};\n" +
" d {};\n" +
" e {};\n" +
" }}\n"
;
std::string rxnTypeFormat = " type {};\n";
std::string rxnEqnFormat = " reaction \"{}\";\n";
#endif