research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Changed the test so that it is more lenient on multiple platforms.

Browse source
This commit is contained in:
Harry Moffat 2008-10-16 01:00:55 +00:00
parent d649d45866
commit 53c4da6ddc
8 changed files with 944 additions and 669 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
test_problems/surfSolverTest/.cvsignore
View file

@ -8,4 +8,10 @@ output2.txt
surfaceSolver
surfaceSolver2
ct2ctml.log
diamond.xml
diff_2csv.txt
diff_2out.txt
diff_csv.txt
diff_out.txt
results.txt
results2.txt

44
test_problems/surfSolverTest/results2_blessed.txt Normal file
View file

@ -0,0 +1,44 @@
Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 H 8.69601e-05 0.00999001 -2.365628e-03
1 N2 0.00773858 0.889011 0.000000e+00
2 CH3 0 0 0.000000e+00
3 CH4 8.69601e-05 0.00999001 0.000000e+00
4 C2H2 8.69601e-05 0.00999001 -1.007796e-04
5 H2 0.00034784 0.03996 1.282503e-03
6 OH 8.69601e-07 9.99001e-05 -7.542496e-05
7 H2O 0.00034784 0.03996 3.880305e-05
8 CO 0 0 3.843405e-05
9 O2 8.69601e-06 0.000999001 -9.060705e-07
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1.4e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293.065 1 1.631252e-04
Bulk Weight Growth Rate = 0.00196 kg/m^2/s
Bulk Growth Rate = 5.57e-07 m/s
Bulk Growth Rate = 2e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0184677 -8.195044e-20
1 Csoot-H 0.981532 -8.195044e-20
Sum of coverages = 1
Surface Phase: soot_interface (11)
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0184677 -8.195044e-20
1 Csoot-H 0.981532 -8.195044e-20
Sum of coverages = 1

148
test_problems/surfSolverTest/results_blessed.txt Normal file
View file

@ -0,0 +1,148 @@
Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 H 8.7e-05 0.00999 -2.366e-03
1 N2 0.00774 0.889 0.000e+00
2 CH3 0 0 0.000e+00
3 CH4 8.7e-05 0.00999 0.000e+00
4 C2H2 8.7e-05 0.00999 -1.008e-04
5 H2 0.000348 0.04 1.283e-03
6 OH 8.7e-07 9.99e-05 -7.542e-05
7 H2O 0.000348 0.04 3.880e-05
8 CO 0 0 3.843e-05
9 O2 8.7e-06 0.000999 -9.061e-07
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1.4e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293 1 1.631e-04
Bulk Weight Growth Rate = 0.00196 kg/m^2/s
Bulk Growth Rate = 5.57e-07 m/s
Bulk Growth Rate = 2e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0185 0.000e+00
1 Csoot-H 0.982 0.000e+00
Sum of coverages = 1
Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 H 0.000113 0.013 -3.122e-03
1 N2 0.00771 0.886 0.000e+00
2 CH3 0 0 0.000e+00
3 CH4 8.7e-05 0.00999 0.000e+00
4 C2H2 8.7e-05 0.00999 -1.015e-04
5 H2 0.000348 0.04 1.661e-03
6 OH 8.7e-07 9.99e-05 -7.541e-05
7 H2O 0.000348 0.04 3.879e-05
8 CO 0 0 3.844e-05
9 O2 8.7e-06 0.000999 -9.129e-07
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1.4e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293 1 1.646e-04
Bulk Weight Growth Rate = 0.00198 kg/m^2/s
Bulk Growth Rate = 5.62e-07 m/s
Bulk Growth Rate = 2.02e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0186 0.000e+00
1 Csoot-H 0.981 0.000e+00
Sum of coverages = 1
Gas Temperature = 1.5e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 H 0.000106 0.013 -3.833e-03
1 N2 0.00722 0.886 0.000e+00
2 CH3 0 0 0.000e+00
3 CH4 8.14e-05 0.00999 0.000e+00
4 C2H2 8.14e-05 0.00999 -1.499e-04
5 H2 0.000326 0.04 2.065e-03
6 OH 8.14e-07 9.99e-05 -7.369e-05
7 H2O 0.000326 0.04 3.846e-05
8 CO 0 0 3.788e-05
9 O2 8.14e-06 0.000999 -1.324e-06
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1.5e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293 1 2.620e-04
Bulk Weight Growth Rate = 0.00315 kg/m^2/s
Bulk Growth Rate = 8.94e-07 m/s
Bulk Growth Rate = 3.22e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1.5e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0243 0.000e+00
1 Csoot-H 0.976 0.000e+00
Sum of coverages = 1
Gas Temperature = 1.5e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 H 0.000106 0.013 -3.833e-03
1 N2 0.00722 0.886 0.000e+00
2 CH3 0 0 0.000e+00
3 CH4 8.14e-05 0.00999 0.000e+00
4 C2H2 8.14e-05 0.00999 -1.499e-04
5 H2 0.000326 0.04 2.065e-03
6 OH 8.14e-07 9.99e-05 -7.369e-05
7 H2O 0.000326 0.04 3.846e-05
8 CO 0 0 3.788e-05
9 O2 8.14e-06 0.000999 -1.324e-06
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1.5e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293 1 2.620e-04
Bulk Weight Growth Rate = 0.00315 kg/m^2/s
Bulk Growth Rate = 8.94e-07 m/s
Bulk Growth Rate = 3.22e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1.5e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0243 0.000e+00
1 Csoot-H 0.976 0.000e+00
Sum of coverages = 1

31
test_problems/surfSolverTest/runtest
View file

@ -3,7 +3,7 @@
#
temp_success="1"
/bin/rm -f output.txt diff_test.txt
/bin/rm -f output.txt diff_test.txt results.txt diff_csv.txt
tname="surfaceSolver"
#################################################################
@ -19,16 +19,31 @@ then
temp_success="0"
echo "surfaceSolver ($tname test) returned with bad status, $retnStat, check output"
fi
diff results.txt results_blessed.txt > diff_csv.txt
retnStat_csv=$?
diff -w output.txt surfaceSolver_blessed.out > diff_test.txt
retnStat=$?
if [ $retnStat = "0" ]
diff -w output.txt surfaceSolver_blessed.out > diff_out.txt
retnStat_txt=$?
if test $retnStat_csv = "0"
then
echo "successful diff comparison on $tname test"
echo "Successful test comparison on "`pwd`
if [ $retnStat_txt != "0" ]
then
echo " But, text output files have differences. See diff_out.txt"
fi
else
echo "unsuccessful diff comparison on $tname test"
echo "FAILED" > csvCode.txt
temp_success="0"
echo "Unsuccessful test comparison on "`pwd` " test"
if test $retnStat_csv != "0"
then
echo " results files are different - see diff_csv.txt"
fi
if test $retnStat_txt != "0"
then
echo " And, text output files have differences. See diff_out.txt"
fi
fi

34
test_problems/surfSolverTest/runtest2
View file

@ -1,10 +1,11 @@
#!/bin/sh
#
#
temp_success="1"
/bin/rm -f output2.txt diff_test2.txt
tname="surfaceSolver2"
#
/bin/rm -f output2.txt diff_2test.txt results2.txt
tname="surfaceSolver"
#################################################################
#
#################################################################
@ -18,16 +19,29 @@ then
temp_success="0"
echo "surfaceSolver ($tname test) returned with bad status, $retnStat, check output"
fi
diff results2.txt results2_blessed.txt > diff_2csv.txt
retnStat_csv=$?
diff -w output2.txt surfaceSolver2_blessed.out > diff_test2.txt
retnStat=$?
if [ $retnStat = "0" ]
diff -w output2.txt surfaceSolver2_blessed.out > diff_2out.txt
retnStat_txt=$?
if test $retnStat_csv = "0"
then
echo "successful diff comparison on $tname test"
echo "Successful test comparison on "`pwd`
if [ $retnStat_txt != "0" ]
then
echo " But, text output files have differences. See diff_2out.txt"
fi
else
echo "unsuccessful diff comparison on $tname test"
echo "FAILED" > csvCode.txt
temp_success="0"
echo "Unsuccessful test comparison on "`pwd` " test"
if test $retnStat_csv != "0"
then
echo " results files are different - see diff_2csv.txt"
fi
if test $retnStat_txt != "0"
then
echo " And, text output files have differences. See diff_2out.txt"
fi
fi

594
test_problems/surfSolverTest/surfaceSolver.cpp
View file

@ -17,6 +17,7 @@
#include <iostream>
#include <string>
#include <vector>
#include <fstream>
#define MSSIZE 200
using namespace std;
@ -43,53 +44,59 @@ static void printUsage()
using namespace Cantera;
void printGas(ThermoPhase *gasTP, InterfaceKinetics * iKin_ptr, double *src) {
void printGas(ostream &oooo, ThermoPhase *gasTP, InterfaceKinetics * iKin_ptr, double *src) {
double x[MSSIZE];
double C[MSSIZE];
char buf[256];
oooo.precision(3);
string gasPhaseName = "gas";
gasTP->getMoleFractions(x);
gasTP->getConcentrations(C);
double Temp = gasTP->temperature();
double p = gasTP->pressure();
cout << "Gas Temperature = " << Temp << endl;
cout << "Gas Pressure = " << p << endl;
oooo << "Gas Temperature = " << Temp << endl;
oooo << "Gas Pressure = " << p << endl;
int kstart = iKin_ptr->kineticsSpeciesIndex(0, 0);
cout << "Gas Phase: " << gasPhaseName << " "
oooo << "Gas Phase: " << gasPhaseName << " "
<< "(" << kstart << ")" << endl;
cout << " Name "
oooo << " Name "
<< " Conc MoleF SrcRate " << endl;
cout << " "
oooo << " "
<< " (kmol/m^3) (kmol/m^2/s) " << endl;
double sum = 0.0;
int nspGas = gasTP->nSpecies();
for (int k = 0; k < nspGas; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(k, 0);
printf("%4d %24s %14g %14g %14e\n",
k, gasTP->speciesName(k).c_str(),
C[k], x[k], src[kstart]);
sprintf(buf, "%4d %24s %14.3g %14.3g %14.3e\n",
k, gasTP->speciesName(k).c_str(),
C[k], x[k], src[kstart]);
oooo << buf;
sum += x[k];
}
cout << "Sum of gas mole fractions= " << sum << endl;
cout << endl;
oooo << "Sum of gas mole fractions= " << sum << endl;
oooo << endl;
}
void printBulk(ThermoPhase *bulkPhaseTP, InterfaceKinetics * iKin_ptr, double *src) {
void printBulk(ostream &oooo,
ThermoPhase *bulkPhaseTP, InterfaceKinetics * iKin_ptr, double *src) {
double x[MSSIZE];
double C[MSSIZE];
char buf[256];
oooo.precision(3);
string bulkParticlePhaseName = bulkPhaseTP->id();
bulkPhaseTP->getMoleFractions(x);
bulkPhaseTP->getConcentrations(C);
int kstart = iKin_ptr->kineticsSpeciesIndex(0, 1);
double dens = bulkPhaseTP->density();
cout << "Bulk Phase: " << bulkParticlePhaseName << " "
oooo << "Bulk Phase: " << bulkParticlePhaseName << " "
<< "(" << kstart << ")" << endl;
double Temp = bulkPhaseTP->temperature();
double p = bulkPhaseTP->pressure();
cout << "Bulk Temperature = " << Temp << endl;
cout << "Bulk Pressure = " << p << endl;
cout << " Name "
oooo << "Bulk Temperature = " << Temp << endl;
oooo << "Bulk Pressure = " << p << endl;
oooo << " Name "
<< " Conc MoleF SrcRate " << endl;
cout << " "
oooo << " "
<< " (kmol/m^3) (kmol/m^2/s) " << endl;
double sum = 0.0;
double Wsum = 0.0;
@ -97,319 +104,342 @@ void printBulk(ThermoPhase *bulkPhaseTP, InterfaceKinetics * iKin_ptr, double *s
int nspBulk = bulkPhaseTP->nSpecies();
for (int k = 0; k < nspBulk; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(k, 1);
printf("%4d %24s %14g %14g %14e\n",
sprintf(buf, "%4d %24s %14.3g %14.3g %14.3e\n",
k, bulkPhaseTP->speciesName(k).c_str(),
C[k], x[k], src[kstart]);
oooo << buf;
sum += x[k];
Wsum += src[kstart] * molecW[k];
}
cout << "Bulk Weight Growth Rate = " << Wsum << " kg/m^2/s" << endl;
oooo.precision(3);
oooo << "Bulk Weight Growth Rate = " << Wsum << " kg/m^2/s" << endl;
double gr = Wsum / dens;
cout << "Bulk Growth Rate = " << gr << " m/s" << endl;
cout << "Bulk Growth Rate = " << gr * 1.0E6 * 3600.
oooo << "Bulk Growth Rate = " << gr << " m/s" << endl;
oooo << "Bulk Growth Rate = " << gr * 1.0E6 * 3600.
<< " microns / hour" << endl;
cout << "Density of bulk phase = " << dens << " kg / m^3 "<< endl;
cout << " = " << dens / 1.0E3
oooo << "Density of bulk phase = " << dens << " kg / m^3 "<< endl;
oooo << " = " << dens / 1.0E3
<<" gm / cm^3 " << endl;
cout << "Sum of bulk mole fractions= " << sum << endl;
cout << endl;
oooo << "Sum of bulk mole fractions= " << sum << endl;
oooo << endl;
}
void printSurf(ThermoPhase *surfPhaseTP, InterfaceKinetics * iKin_ptr, double *src) {
void printSurf(ostream &oooo,
ThermoPhase *surfPhaseTP, InterfaceKinetics * iKin_ptr, double *src) {
double x[MSSIZE];
char buf[256];
string surfParticlePhaseName = surfPhaseTP->id();
surfPhaseTP->getMoleFractions(x);
int kstart = iKin_ptr->kineticsSpeciesIndex(0, 2);
cout << "Surface Phase: " << surfParticlePhaseName
oooo << "Surface Phase: " << surfParticlePhaseName
<< " (" << kstart << ")" << endl;
double Temp = surfPhaseTP->temperature();
double p = surfPhaseTP->pressure();
cout << "Surface Temperature = " << Temp << endl;
cout << "Surface Pressure = " << p << endl;
cout << " Name "
oooo << "Surface Temperature = " << Temp << endl;
oooo << "Surface Pressure = " << p << endl;
oooo << " Name "
<< " Coverage SrcRate " << endl;
double sum = 0.0;
int nspSurf = surfPhaseTP->nSpecies();
int nspSurf = surfPhaseTP->nSpecies();
for (int k = 0; k < nspSurf; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(0, 2);
printf("%4d %24s %14g %14e\n",
double srcK = src[kstart];
if (fabs(srcK) < 1.0E-7) {
srcK = 0.0;
}
sprintf(buf, "%4d %24s %14.3g %14.3e\n",
k, surfPhaseTP->speciesName(k).c_str(),
x[k], src[kstart]);
x[k], srcK);
oooo << buf;
sum += x[k];
}
cout << "Sum of coverages = " << sum << endl;
oooo << "Sum of coverages = " << sum << endl;
}
int main(int argc, char** argv) {
string infile;
int ioflag = 1;
int i, k;
// look for command-line options
if (argc > 1) {
string tok;
for (int j = 1; j < argc; j++) {
tok = string(argv[j]);
if (tok[0] == '-') {
int nopt = static_cast<int>(tok.size());
for (int n = 1; n < nopt; n++) {
if (tok[n] == 'h') {
printUsage();
exit(0);
} else if (tok[n] == 'd') {
int lvl = 0;
if (j < (argc - 1)) {
string tokla = string(argv[j+1]);
if (strlen(tokla.c_str()) > 0) {
lvl = atoi(tokla.c_str());
n = nopt - 1;
j += 1;
ioflag = lvl;
}
}
} else {
printUsage();
exit(1);
string infile;
int ioflag = 1;
int i, k;
// look for command-line options
if (argc > 1) {
string tok;
for (int j = 1; j < argc; j++) {
tok = string(argv[j]);
if (tok[0] == '-') {
int nopt = static_cast<int>(tok.size());
for (int n = 1; n < nopt; n++) {
if (tok[n] == 'h') {
printUsage();
exit(0);
} else if (tok[n] == 'd') {
int lvl = 0;
if (j < (argc - 1)) {
string tokla = string(argv[j+1]);
if (strlen(tokla.c_str()) > 0) {
lvl = atoi(tokla.c_str());
n = nopt - 1;
j += 1;
ioflag = lvl;
}
}
}
} else if (infile == "") {
infile = tok;
}
else {
printUsage();
exit(1);
}
}
}
if (infile == "") {
infile = "diamond.cti";
}
try {
/*************************************************************/
/*
* FILL IN THESE NAMES FOR EACH PROBLEM
*/
/*
* ProblemNumber = 0 : diamond.cti
* = 1 : haca.cti
*/
int ProblemNumber = 1;
string gasPhaseName = "gas";
string bulkParticlePhaseName = "diamond";
string surfParticlePhaseName = "diamond_100";
if (ProblemNumber == 1) {
gasPhaseName = "gas";
bulkParticlePhaseName = "soot";
surfParticlePhaseName = "soot_interface";
}
/************************************************************/
XML_Node *xc = new XML_Node();
string path = findInputFile(infile);
ctml::get_CTML_Tree(xc, path);
XML_Node * const xg = (XML_Node *) findXMLPhase(xc, gasPhaseName);
if (!xg) {
printf("ERROR: Could not find gas phase named, %s, in file\n",
gasPhaseName.c_str());
exit(-1);
}
ThermoPhase *gasTP = newPhase(*xg);
int nspGas = gasTP->nSpecies();
cout << "Number of species = " << nspGas << endl;
XML_Node * const xd =
(XML_Node *) findXMLPhase(xc, bulkParticlePhaseName);
if (!xd) {
printf("ERROR: Could not find bulk phase named, %s, in file\n",
bulkParticlePhaseName.c_str());
exit(-1);
}
ThermoPhase *bulkPhaseTP = newPhase(*xd);
int nspBulk = bulkPhaseTP->nSpecies();
cout << "Number of species in bulk phase named " <<
bulkParticlePhaseName << " = " << nspBulk << endl;
XML_Node * const xs =
(XML_Node *) findXMLPhase(xc, surfParticlePhaseName);
if (!xs) {
printf("ERROR: Could not find surf Particle phase named, %s, in file\n",
surfParticlePhaseName.c_str());
exit(-1);
}
ThermoPhase *surfPhaseTP = newPhase(*xs);
int nsp_d100 = surfPhaseTP->nSpecies();
cout << "Number of species in surface phase, " << surfParticlePhaseName
<< " = " << nsp_d100 << endl;
vector<ThermoPhase *> phaseList;
phaseList.push_back(gasTP);
phaseList.push_back(bulkPhaseTP);
phaseList.push_back(surfPhaseTP);
InterfaceKinetics *iKin_ptr = new InterfaceKinetics();
importKinetics(*xs, phaseList, iKin_ptr);
int nr = iKin_ptr->nReactions();
cout << "Number of reactions = " << nr << endl;
double x[MSSIZE], p = OneAtm;
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
x[0] = 0.0010;
x[1] = 0.9888;
x[2] = 0.0002;
x[3] = 0.0100;
p = 20.0*OneAtm/760.0;
gasTP->setState_TPX(1200., p, x);
}
/*
* Set the surface initial state
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
int i0 = surfPhaseTP->speciesIndex("c6H*");
if (i0 >= 0) {
x[i0] = 0.1;
}
int i1 = surfPhaseTP->speciesIndex("c6HH");
if (i1 >= 0) {
x[i1] = 0.9;
}
surfPhaseTP->setState_TX(1200., x);
}
/*
* Set the bulk Phase State
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
x[0] = 1.0;
bulkPhaseTP->setState_TPX(1200., p, x);
}
iKin_ptr->setIOFlag(ioflag);
/*
* Solve the Equation system
*/
//iKin_ptr->advanceCoverages(100.);
iKin_ptr->solvePseudoSteadyStateProblem();
/*
* Download the source terms for the rate equations
*/
double src[MSSIZE];
iKin_ptr->getNetProductionRates(src);
double sum = 0.0;
if (ProblemNumber == 0) {
double naH;
for (k = 0; k < 13; k++) {
if (k < 4) {
naH = gasTP->nAtoms(k, 0);
} else if (k == 4) {
naH = 0;
} else if (k > 4) {
int itp = k - 5;
naH = surfPhaseTP->nAtoms(itp, 0);
}
cout << k << " " << naH << " " ;
if (fabs(src[k]) < 2.0E-17) {
cout << " nil" << endl;
} else {
cout << src[k] << endl;
printUsage();
exit(1);
}
sum += naH * src[k];
}
cout << "sum = " << sum << endl;
} else if (infile == "") {
infile = tok;
}
else {
printUsage();
exit(1);
}
}
}
if (infile == "") {
infile = "diamond.cti";
}
try {
/*************************************************************/
/*
* FILL IN THESE NAMES FOR EACH PROBLEM
*/
/*
* ProblemNumber = 0 : diamond.cti
* = 1 : haca.cti
*/
int ProblemNumber = 1;
string gasPhaseName = "gas";
string bulkParticlePhaseName = "diamond";
string surfParticlePhaseName = "diamond_100";
if (ProblemNumber == 1) {
gasPhaseName = "gas";
bulkParticlePhaseName = "soot";
surfParticlePhaseName = "soot_interface";
}
/************************************************************/
XML_Node *xc = new XML_Node();
string path = findInputFile(infile);
ctml::get_CTML_Tree(xc, path);
XML_Node * const xg = (XML_Node *) findXMLPhase(xc, gasPhaseName);
if (!xg) {
printf("ERROR: Could not find gas phase named, %s, in file\n",
gasPhaseName.c_str());
exit(-1);
}
ThermoPhase *gasTP = newPhase(*xg);
int nspGas = gasTP->nSpecies();
cout << "Number of species = " << nspGas << endl;
XML_Node * const xd =
(XML_Node *) findXMLPhase(xc, bulkParticlePhaseName);
if (!xd) {
printf("ERROR: Could not find bulk phase named, %s, in file\n",
bulkParticlePhaseName.c_str());
exit(-1);
}
ThermoPhase *bulkPhaseTP = newPhase(*xd);
int nspBulk = bulkPhaseTP->nSpecies();
cout << "Number of species in bulk phase named " <<
bulkParticlePhaseName << " = " << nspBulk << endl;
XML_Node * const xs =
(XML_Node *) findXMLPhase(xc, surfParticlePhaseName);
if (!xs) {
printf("ERROR: Could not find surf Particle phase named, %s, in file\n",
surfParticlePhaseName.c_str());
exit(-1);
}
ThermoPhase *surfPhaseTP = newPhase(*xs);
int nsp_d100 = surfPhaseTP->nSpecies();
cout << "Number of species in surface phase, " << surfParticlePhaseName
<< " = " << nsp_d100 << endl;
vector<ThermoPhase *> phaseList;
phaseList.push_back(gasTP);
phaseList.push_back(bulkPhaseTP);
phaseList.push_back(surfPhaseTP);
InterfaceKinetics *iKin_ptr = new InterfaceKinetics();
importKinetics(*xs, phaseList, iKin_ptr);
int nr = iKin_ptr->nReactions();
cout << "Number of reactions = " << nr << endl;
double x[MSSIZE], p = OneAtm;
ofstream ofile("results.txt");
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
x[0] = 0.0010;
x[1] = 0.9888;
x[2] = 0.0002;
x[3] = 0.0100;
p = 20.0*OneAtm/760.0;
gasTP->setState_TPX(1200., p, x);
}
/*
* Set the surface initial state
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
int i0 = surfPhaseTP->speciesIndex("c6H*");
if (i0 >= 0) {
x[i0] = 0.1;
}
int i1 = surfPhaseTP->speciesIndex("c6HH");
if (i1 >= 0) {
x[i1] = 0.9;
}
surfPhaseTP->setState_TX(1200., x);
}
/*
* Set the bulk Phase State
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
x[0] = 1.0;
bulkPhaseTP->setState_TPX(1200., p, x);
}
iKin_ptr->setIOFlag(ioflag);
/*
* Solve the Equation system
*/
//iKin_ptr->advanceCoverages(100.);
iKin_ptr->solvePseudoSteadyStateProblem();
/*
* Download the source terms for the rate equations
*/
double src[MSSIZE];
iKin_ptr->getNetProductionRates(src);
double sum = 0.0;
if (ProblemNumber == 0) {
double naH;
for (k = 0; k < 13; k++) {
if (k < 4) {
naH = gasTP->nAtoms(k, 0);
} else if (k == 4) {
naH = 0;
} else if (k > 4) {
int itp = k - 5;
naH = surfPhaseTP->nAtoms(itp, 0);
}
cout << k << " " << naH << " " ;
if (fabs(src[k]) < 2.0E-17) {
cout << " nil" << endl;
} else {
cout << src[k] << endl;
}
sum += naH * src[k];
}
cout << "sum = " << sum << endl;
}
printGas(gasTP, iKin_ptr, src);
printBulk(bulkPhaseTP, iKin_ptr, src);
printSurf(surfPhaseTP, iKin_ptr, src) ;
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
printGas(ofile, gasTP, iKin_ptr, src);
printBulk(ofile, bulkPhaseTP, iKin_ptr, src);
printSurf(ofile, surfPhaseTP, iKin_ptr, src) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
double pres = gasTP->pressure();
gasTP->getMoleFractions(x);
double tmp = 0.3 * x[0];
double tmp2 = 0.3 * x[1];
if (tmp2 < tmp) tmp = tmp2;
x[0] += tmp;
x[1] -= tmp;
gasTP->setState_PX(pres, x);
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
double pres = gasTP->pressure();
gasTP->getMoleFractions(x);
double tmp = 0.3 * x[0];
double tmp2 = 0.3 * x[1];
if (tmp2 < tmp) tmp = tmp2;
x[0] += tmp;
x[1] -= tmp;
gasTP->setState_PX(pres, x);
iKin_ptr->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin_ptr->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
printGas(gasTP, iKin_ptr, src);
printBulk(bulkPhaseTP, iKin_ptr, src);
printSurf(surfPhaseTP, iKin_ptr, src) ;
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printGas(ofile, gasTP, iKin_ptr, src);
printBulk(ofile, bulkPhaseTP, iKin_ptr, src);
printSurf(ofile, surfPhaseTP, iKin_ptr, src) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
pres = gasTP->pressure();
double temp = gasTP->temperature();
temp += 95;
gasTP->setState_TP(temp, pres);
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
pres = gasTP->pressure();
double temp = gasTP->temperature();
temp += 95;
gasTP->setState_TP(temp, pres);
iKin_ptr->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin_ptr->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
printGas(gasTP, iKin_ptr, src);
printBulk(bulkPhaseTP, iKin_ptr, src);
printSurf(surfPhaseTP, iKin_ptr, src) ;
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printGas(ofile, gasTP, iKin_ptr, src);
printBulk(ofile, bulkPhaseTP, iKin_ptr, src);
printSurf(ofile, surfPhaseTP, iKin_ptr, src) ;
/*****************************************************************************/
/* Now Don't Tweak the inputs at all */
/****************************************************************************/
gasTP->setState_TP(temp, pres);
/*****************************************************************************/
/* Now Don't Tweak the inputs at all */
/****************************************************************************/
gasTP->setState_TP(temp, pres);
iKin_ptr->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin_ptr->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
printGas(gasTP, iKin_ptr, src);
printBulk(bulkPhaseTP, iKin_ptr, src);
printSurf(surfPhaseTP, iKin_ptr, src) ;
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
delete iKin_ptr;
delete gasTP; gasTP = 0;
delete bulkPhaseTP; bulkPhaseTP = 0;
delete surfPhaseTP; surfPhaseTP = 0;
delete xc;
appdelete();
}
catch (CanteraError) {
showErrors(cout);
}
printGas(ofile, gasTP, iKin_ptr, src);
printBulk(ofile, bulkPhaseTP, iKin_ptr, src);
printSurf(ofile, surfPhaseTP, iKin_ptr, src) ;
return 0;
delete iKin_ptr;
delete gasTP; gasTP = 0;
delete bulkPhaseTP; bulkPhaseTP = 0;
delete surfPhaseTP; surfPhaseTP = 0;
delete xc;
appdelete();
}
catch (CanteraError) {
showErrors(cout);
}
return 0;
}
/***********************************************************/

658
test_problems/surfSolverTest/surfaceSolver2.cpp
View file

@ -45,55 +45,60 @@ using namespace Cantera;
void printGas(ThermoPhase *gasTP, InterfaceKinetics * iKin_ptr, double *src) {
void printGas(ostream &oooo, ThermoPhase *gasTP, InterfaceKinetics * iKin_ptr, double *src) {
double x[MSSIZE];
double C[MSSIZE];
char buf[256];
oooo.precision(3);
string gasPhaseName = "gas";
gasTP->getMoleFractions(x);
gasTP->getConcentrations(C);
double Temp = gasTP->temperature();
double p = gasTP->pressure();
cout << "Gas Temperature = " << Temp << endl;
cout << "Gas Pressure = " << p << endl;
oooo << "Gas Temperature = " << Temp << endl;
oooo << "Gas Pressure = " << p << endl;
int kstart = iKin_ptr->kineticsSpeciesIndex(0, 0);
cout << "Gas Phase: " << gasPhaseName << " "
oooo << "Gas Phase: " << gasPhaseName << " "
<< "(" << kstart << ")" << endl;
cout << " Name "
oooo << " Name "
<< " Conc MoleF SrcRate " << endl;
cout << " "
oooo << " "
<< " (kmol/m^3) (kmol/m^2/s) " << endl;
double sum = 0.0;
int nspGas = gasTP->nSpecies();
for (int k = 0; k < nspGas; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(k, 0);
printf("%4d %24s %14g %14g %14e\n",
k, gasTP->speciesName(k).c_str(),
C[k], x[k], src[kstart]);
sprintf(buf, "%4d %24s %14g %14g %14e\n",
k, gasTP->speciesName(k).c_str(),
C[k], x[k], src[kstart]);
oooo << buf;
sum += x[k];
}
cout << "Sum of gas mole fractions= " << sum << endl;
cout << endl;
oooo << "Sum of gas mole fractions= " << sum << endl;
oooo << endl;
}
void printBulk(ThermoPhase *bulkPhaseTP, InterfaceKinetics * iKin_ptr, double *src) {
void printBulk(ostream &oooo,ThermoPhase *bulkPhaseTP, InterfaceKinetics * iKin_ptr, double *src) {
double x[MSSIZE];
double C[MSSIZE];
char buf[256];
oooo.precision(3);
string bulkParticlePhaseName = bulkPhaseTP->id();
bulkPhaseTP->getMoleFractions(x);
bulkPhaseTP->getConcentrations(C);
int kstart = iKin_ptr->kineticsSpeciesIndex(0, 1);
double dens = bulkPhaseTP->density();
cout << "Bulk Phase: " << bulkParticlePhaseName << " "
oooo << "Bulk Phase: " << bulkParticlePhaseName << " "
<< "(" << kstart << ")" << endl;
double Temp = bulkPhaseTP->temperature();
double p = bulkPhaseTP->pressure();
cout << "Bulk Temperature = " << Temp << endl;
cout << "Bulk Pressure = " << p << endl;
cout << " Name "
oooo << "Bulk Temperature = " << Temp << endl;
oooo << "Bulk Pressure = " << p << endl;
oooo << " Name "
<< " Conc MoleF SrcRate " << endl;
cout << " "
oooo << " "
<< " (kmol/m^3) (kmol/m^2/s) " << endl;
double sum = 0.0;
double Wsum = 0.0;
@ -101,362 +106,375 @@ void printBulk(ThermoPhase *bulkPhaseTP, InterfaceKinetics * iKin_ptr, double *s
int nspBulk = bulkPhaseTP->nSpecies();
for (int k = 0; k < nspBulk; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(k, 1);
printf("%4d %24s %14g %14g %14e\n",
sprintf(buf, "%4d %24s %14g %14g %14e\n",
k, bulkPhaseTP->speciesName(k).c_str(),
C[k], x[k], src[kstart]);
oooo << buf;
sum += x[k];
Wsum += src[kstart] * molecW[k];
}
cout << "Bulk Weight Growth Rate = " << Wsum << " kg/m^2/s" << endl;
oooo << "Bulk Weight Growth Rate = " << Wsum << " kg/m^2/s" << endl;
double gr = Wsum / dens;
cout << "Bulk Growth Rate = " << gr << " m/s" << endl;
cout << "Bulk Growth Rate = " << gr * 1.0E6 * 3600.
oooo << "Bulk Growth Rate = " << gr << " m/s" << endl;
oooo << "Bulk Growth Rate = " << gr * 1.0E6 * 3600.
<< " microns / hour" << endl;
cout << "Density of bulk phase = " << dens << " kg / m^3 "<< endl;
cout << " = " << dens / 1.0E3
oooo << "Density of bulk phase = " << dens << " kg / m^3 "<< endl;
oooo << " = " << dens / 1.0E3
<<" gm / cm^3 " << endl;
cout << "Sum of bulk mole fractions= " << sum << endl;
cout << endl;
oooo << "Sum of bulk mole fractions= " << sum << endl;
oooo << endl;
}
void printSurf(ThermoPhase *surfPhaseTP, InterfaceKinetics * iKin_ptr, double *src) {
void printSurf(ostream &oooo, ThermoPhase *surfPhaseTP,
InterfaceKinetics * iKin_ptr, double *src) {
double x[MSSIZE];
char buf[256];
oooo.precision(3);
string surfParticlePhaseName = surfPhaseTP->id();
surfPhaseTP->getMoleFractions(x);
int kstart = iKin_ptr->kineticsSpeciesIndex(0, 2);
cout << "Surface Phase: " << surfParticlePhaseName
oooo << "Surface Phase: " << surfParticlePhaseName
<< " (" << kstart << ")" << endl;
double Temp = surfPhaseTP->temperature();
double p = surfPhaseTP->pressure();
cout << "Surface Temperature = " << Temp << endl;
cout << "Surface Pressure = " << p << endl;
cout << " Name "
oooo << "Surface Temperature = " << Temp << endl;
oooo << "Surface Pressure = " << p << endl;
oooo << " Name "
<< " Coverage SrcRate " << endl;
double sum = 0.0;
int nspSurf = surfPhaseTP->nSpecies();
for (int k = 0; k < nspSurf; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(0, 2);
printf("%4d %24s %14g %14e\n",
sprintf(buf, "%4d %24s %14g %14e\n",
k, surfPhaseTP->speciesName(k).c_str(),
x[k], src[kstart]);
oooo << buf;
sum += x[k];
}
cout << "Sum of coverages = " << sum << endl;
oooo << "Sum of coverages = " << sum << endl;
}
int main(int argc, char** argv) {
string infile;
int i, k;
int ioflag = 1;
// look for command-line options
if (argc > 1) {
string tok;
for (int j = 1; j < argc; j++) {
tok = string(argv[j]);
if (tok[0] == '-') {
int nopt = static_cast<int>(tok.size());
for (int n = 1; n < nopt; n++) {
if (tok[n] == 'h') {
printUsage();
exit(0);
} else if (tok[n] == 'd') {
int lvl = 0;
if (j < (argc - 1)) {
string tokla = string(argv[j+1]);
if (strlen(tokla.c_str()) > 0) {
lvl = atoi(tokla.c_str());
n = nopt - 1;
j += 1;
ioflag = lvl;
}
string infile;
int i, k;
int ioflag = 1;
// look for command-line options
if (argc > 1) {
string tok;
for (int j = 1; j < argc; j++) {
tok = string(argv[j]);
if (tok[0] == '-') {
int nopt = static_cast<int>(tok.size());
for (int n = 1; n < nopt; n++) {
if (tok[n] == 'h') {
printUsage();
exit(0);
} else if (tok[n] == 'd') {
int lvl = 0;
if (j < (argc - 1)) {
string tokla = string(argv[j+1]);
if (strlen(tokla.c_str()) > 0) {
lvl = atoi(tokla.c_str());
n = nopt - 1;
j += 1;
ioflag = lvl;
}
} else {
printUsage();
exit(1);
}
} else {
printUsage();
exit(1);
}
} else if (infile == "") {
infile = tok;
}
else {
printUsage();
exit(1);
}
} else if (infile == "") {
infile = tok;
}
else {
printUsage();
exit(1);
}
}
if (infile == "") {
infile = "diamond.cti";
}
}
if (infile == "") {
infile = "diamond.cti";
}
try {
/*************************************************************/
try {
/*************************************************************/
/*
* FILL IN THESE NAMES FOR EACH PROBLEM
*/
/*
* ProblemNumber = 0 : diamond.cti
* = 1 : haca.cti
*/
int ProblemNumber = 1;
string gasPhaseName = "gas";
string bulkParticlePhaseName = "diamond";
string surfParticlePhaseName = "diamond_100";
if (ProblemNumber == 1) {
gasPhaseName = "gas";
bulkParticlePhaseName = "soot";
surfParticlePhaseName = "soot_interface";
}
/*
* FILL IN THESE NAMES FOR EACH PROBLEM
*/
/*
* ProblemNumber = 0 : diamond.cti
* = 1 : haca.cti
*/
int ProblemNumber = 1;
string gasPhaseName = "gas";
string bulkParticlePhaseName = "diamond";
string surfParticlePhaseName = "diamond_100";
if (ProblemNumber == 1) {
gasPhaseName = "gas";
bulkParticlePhaseName = "soot";
surfParticlePhaseName = "soot_interface";
}
/************************************************************/
XML_Node *xc = new XML_Node();
string path = findInputFile(infile);
ctml::get_CTML_Tree(xc, path);
/************************************************************/
XML_Node *xc = new XML_Node();
string path = findInputFile(infile);
ctml::get_CTML_Tree(xc, path);
XML_Node * const xg = (XML_Node *) findXMLPhase(xc, gasPhaseName);
if (!xg) {
printf("ERROR: Could not find gas phase named, %s, in file\n",
gasPhaseName.c_str());
exit(-1);
}
ThermoPhase *gasTP = newPhase(*xg);
int nspGas = gasTP->nSpecies();
cout << "Number of species = " << nspGas << endl;
XML_Node * const xg = (XML_Node *) findXMLPhase(xc, gasPhaseName);
if (!xg) {
printf("ERROR: Could not find gas phase named, %s, in file\n",
gasPhaseName.c_str());
exit(-1);
}
ThermoPhase *gasTP = newPhase(*xg);
int nspGas = gasTP->nSpecies();
cout << "Number of species = " << nspGas << endl;
XML_Node * const xd =
(XML_Node *) findXMLPhase(xc, bulkParticlePhaseName);
if (!xd) {
printf("ERROR: Could not find bulk phase named, %s, in file\n",
bulkParticlePhaseName.c_str());
exit(-1);
}
ThermoPhase *bulkPhaseTP = newPhase(*xd);
int nspBulk = bulkPhaseTP->nSpecies();
cout << "Number of species in bulk phase named " <<
bulkParticlePhaseName << " = " << nspBulk << endl;
XML_Node * const xd =
(XML_Node *) findXMLPhase(xc, bulkParticlePhaseName);
if (!xd) {
printf("ERROR: Could not find bulk phase named, %s, in file\n",
bulkParticlePhaseName.c_str());
exit(-1);
}
ThermoPhase *bulkPhaseTP = newPhase(*xd);
int nspBulk = bulkPhaseTP->nSpecies();
cout << "Number of species in bulk phase named " <<
bulkParticlePhaseName << " = " << nspBulk << endl;
XML_Node * const xs =
(XML_Node *) findXMLPhase(xc, surfParticlePhaseName);
if (!xs) {
printf("ERROR: Could not find surf Particle phase named,"
"%s, in file\n",
surfParticlePhaseName.c_str());
exit(-1);
}
ThermoPhase *surfPhaseTP = newPhase(*xs);
int nsp_d100 = surfPhaseTP->nSpecies();
cout << "Number of species in surface phase, " << surfParticlePhaseName
<< " = " << nsp_d100 << endl;
XML_Node * const xs =
(XML_Node *) findXMLPhase(xc, surfParticlePhaseName);
if (!xs) {
printf("ERROR: Could not find surf Particle phase named,"
"%s, in file\n",
surfParticlePhaseName.c_str());
exit(-1);
}
ThermoPhase *surfPhaseTP = newPhase(*xs);
int nsp_d100 = surfPhaseTP->nSpecies();
cout << "Number of species in surface phase, " << surfParticlePhaseName
<< " = " << nsp_d100 << endl;
vector<ThermoPhase *> phaseList;
phaseList.push_back(gasTP);
phaseList.push_back(bulkPhaseTP);
phaseList.push_back(surfPhaseTP);
vector<ThermoPhase *> phaseList;
phaseList.push_back(gasTP);
phaseList.push_back(bulkPhaseTP);
phaseList.push_back(surfPhaseTP);
InterfaceKinetics *iKin_ptr = new InterfaceKinetics();
importKinetics(*xs, phaseList, iKin_ptr);
int nr = iKin_ptr->nReactions();
cout << "Number of reactions = " << nr << endl;
InterfaceKinetics *iKin_ptr = new InterfaceKinetics();
importKinetics(*xs, phaseList, iKin_ptr);
int nr = iKin_ptr->nReactions();
cout << "Number of reactions = " << nr << endl;
ofstream ofile("results2.txt");
// create a second copy of the same surface phase
// (this is a made up problem btw to check the software capability)
ThermoPhase *surfPhaseTP2 = newPhase(*xs);
int nsp2 = surfPhaseTP2->nSpecies();
string pname = surfPhaseTP2->id();
cout << "Number of species in 2nd surface phase, " << pname
<< " = " << nsp2 << endl;
// create a second copy of the same surface phase
// (this is a made up problem btw to check the software capability)
ThermoPhase *surfPhaseTP2 = newPhase(*xs);
int nsp2 = surfPhaseTP2->nSpecies();
string pname = surfPhaseTP2->id();
cout << "Number of species in 2nd surface phase, " << pname
<< " = " << nsp2 << endl;
vector<ThermoPhase *> phaseList2;
phaseList2.push_back(gasTP);
phaseList2.push_back(bulkPhaseTP);
phaseList2.push_back(surfPhaseTP2);
vector<ThermoPhase *> phaseList2;
phaseList2.push_back(gasTP);
phaseList2.push_back(bulkPhaseTP);
phaseList2.push_back(surfPhaseTP2);
// create the second InterfaceKinetics object based on the
// second surface phase.
InterfaceKinetics *iKin2_ptr = new InterfaceKinetics();
importKinetics(*xs, phaseList2, iKin2_ptr);
nr = iKin_ptr->nReactions();
cout << "Number of reactions = " << nr << endl;
// create the second InterfaceKinetics object based on the
// second surface phase.
InterfaceKinetics *iKin2_ptr = new InterfaceKinetics();
importKinetics(*xs, phaseList2, iKin2_ptr);
nr = iKin_ptr->nReactions();
cout << "Number of reactions = " << nr << endl;
double x[MSSIZE], p = OneAtm;
double x[MSSIZE], p = OneAtm;
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
x[0] = 0.0010;
x[1] = 0.9888;
x[2] = 0.0002;
x[3] = 0.0100;
p = 20.0*OneAtm/760.0;
gasTP->setState_TPX(1200., p, x);
}
/*
* Set the surface initial state
* other problem numbers take their initial state from the xml files.
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
int i0 = surfPhaseTP->speciesIndex("c6H*");
if (i0 >= 0) {
x[i0] = 0.1;
}
int i1 = surfPhaseTP->speciesIndex("c6HH");
if (i1 >= 0) {
x[i1] = 0.9;
}
surfPhaseTP->setState_TX(1200., x);
}
/*
* Set the bulk Phase State
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
x[0] = 1.0;
bulkPhaseTP->setState_TPX(1200., p, x);
}
/*
* Set-up the Surface Problem
* This problem will consist of 2 identical InterfaceKinetics objects
*/
vector<InterfaceKinetics*> vecKinPtrs;
vecKinPtrs.push_back(iKin_ptr);
vecKinPtrs.push_back(iKin2_ptr);
// Create the ImplicitSurfChem problem
// Initialize it and call the pseudo steadystate capability.
ImplicitSurfChem *surfaceProb = new ImplicitSurfChem(vecKinPtrs);
surfaceProb->initialize();
surfaceProb->setIOFlag(ioflag);
surfaceProb->solvePseudoSteadyStateProblem();
/*
* Download the source terms for the rate equations
*/
double src[MSSIZE];
double src2[MSSIZE];
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
double sum = 0.0;
if (ProblemNumber == 0) {
double naH;
for (k = 0; k < 13; k++) {
if (k < 4) {
naH = gasTP->nAtoms(k, 0);
} else if (k == 4) {
naH = 0;
} else if (k > 4) {
int itp = k - 5;
naH = surfPhaseTP->nAtoms(itp, 0);
}
cout << k << " " << naH << " " ;
if (fabs(src[k]) < 2.0E-17) {
cout << " nil" << endl;
} else {
cout << src[k] << endl;
}
sum += naH * src[k];
}
cout << "sum = " << sum << endl;
}
printGas(gasTP, iKin_ptr, src);
printBulk(bulkPhaseTP, iKin_ptr, src);
printSurf(surfPhaseTP, iKin_ptr, src) ;
printSurf(surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
double pres = gasTP->pressure();
gasTP->getMoleFractions(x);
double tmp = 0.3 * x[0];
double tmp2 = 0.3 * x[1];
if (tmp2 < tmp) tmp = tmp2;
x[0] += tmp;
x[1] -= tmp;
gasTP->setState_PX(pres, x);
surfaceProb->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(gasTP, iKin_ptr, src);
printBulk(bulkPhaseTP, iKin_ptr, src);
printSurf(surfPhaseTP, iKin_ptr, src) ;
printSurf(surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
pres = gasTP->pressure();
double temp = gasTP->temperature();
temp += 95;
gasTP->setState_TP(temp, pres);
surfaceProb->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(gasTP, iKin_ptr, src);
printBulk(bulkPhaseTP, iKin_ptr, src);
printSurf(surfPhaseTP, iKin_ptr, src) ;
printSurf(surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Don't Tweak the inputs at all */
/****************************************************************************/
gasTP->setState_TP(temp, pres);
surfaceProb->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(gasTP, iKin_ptr, src);
printBulk(bulkPhaseTP, iKin_ptr, src);
printSurf(surfPhaseTP, iKin_ptr, src) ;
printSurf(surfPhaseTP2, iKin2_ptr, src2) ;
delete surfaceProb; surfaceProb = 0;
delete iKin_ptr; iKin_ptr = 0;
delete iKin2_ptr; iKin2_ptr = 0;
delete gasTP; gasTP = 0;
delete bulkPhaseTP; bulkPhaseTP = 0;
delete surfPhaseTP; surfPhaseTP = 0;
delete surfPhaseTP2; surfPhaseTP2 = 0;
delete xc;
appdelete();
}
catch (CanteraError) {
showErrors(cout);
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
x[0] = 0.0010;
x[1] = 0.9888;
x[2] = 0.0002;
x[3] = 0.0100;
p = 20.0*OneAtm/760.0;
gasTP->setState_TPX(1200., p, x);
}
return 0;
/*
* Set the surface initial state
* other problem numbers take their initial state from the xml files.
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
int i0 = surfPhaseTP->speciesIndex("c6H*");
if (i0 >= 0) {
x[i0] = 0.1;
}
int i1 = surfPhaseTP->speciesIndex("c6HH");
if (i1 >= 0) {
x[i1] = 0.9;
}
surfPhaseTP->setState_TX(1200., x);
}
/*
* Set the bulk Phase State
*/
for (i = 0; i < MSSIZE; i++) x[i] = 0.0;
if (ProblemNumber == 0) {
x[0] = 1.0;
bulkPhaseTP->setState_TPX(1200., p, x);
}
/*
* Set-up the Surface Problem
* This problem will consist of 2 identical InterfaceKinetics objects
*/
vector<InterfaceKinetics*> vecKinPtrs;
vecKinPtrs.push_back(iKin_ptr);
vecKinPtrs.push_back(iKin2_ptr);
// Create the ImplicitSurfChem problem
// Initialize it and call the pseudo steadystate capability.
ImplicitSurfChem *surfaceProb = new ImplicitSurfChem(vecKinPtrs);
surfaceProb->initialize();
surfaceProb->setIOFlag(ioflag);
surfaceProb->solvePseudoSteadyStateProblem();
/*
* Download the source terms for the rate equations
*/
double src[MSSIZE];
double src2[MSSIZE];
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
double sum = 0.0;
if (ProblemNumber == 0) {
double naH;
for (k = 0; k < 13; k++) {
if (k < 4) {
naH = gasTP->nAtoms(k, 0);
} else if (k == 4) {
naH = 0;
} else if (k > 4) {
int itp = k - 5;
naH = surfPhaseTP->nAtoms(itp, 0);
}
cout << k << " " << naH << " " ;
if (fabs(src[k]) < 2.0E-17) {
cout << " nil" << endl;
} else {
cout << src[k] << endl;
}
sum += naH * src[k];
}
cout << "sum = " << sum << endl;
}
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ;
printGas(ofile, gasTP, iKin_ptr, src);
printBulk(ofile, bulkPhaseTP, iKin_ptr, src);
printSurf(ofile, surfPhaseTP, iKin_ptr, src) ;
printSurf(ofile, surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
double pres = gasTP->pressure();
gasTP->getMoleFractions(x);
double tmp = 0.3 * x[0];
double tmp2 = 0.3 * x[1];
if (tmp2 < tmp) tmp = tmp2;
x[0] += tmp;
x[1] -= tmp;
gasTP->setState_PX(pres, x);
surfaceProb->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
/*
* Set the Gas State:
* -> note that the states are set in the xml files too
*/
pres = gasTP->pressure();
double temp = gasTP->temperature();
temp += 95;
gasTP->setState_TP(temp, pres);
surfaceProb->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Don't Tweak the inputs at all */
/****************************************************************************/
gasTP->setState_TP(temp, pres);
surfaceProb->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ;
delete surfaceProb; surfaceProb = 0;
delete iKin_ptr; iKin_ptr = 0;
delete iKin2_ptr; iKin2_ptr = 0;
delete gasTP; gasTP = 0;
delete bulkPhaseTP; bulkPhaseTP = 0;
delete surfPhaseTP; surfPhaseTP = 0;
delete surfPhaseTP2; surfPhaseTP2 = 0;
delete xc;
appdelete();
}
catch (CanteraError) {
showErrors(cout);
}
return 0;
}
/***********************************************************/

96
test_problems/surfSolverTest/surfaceSolver2_blessed.out
View file

@ -24,8 +24,8 @@ Number of reactions = 8
6 4.7372e+00 4.7359e+00 5.8161e-01 2.9647e-01 Csoot-*
7 6.1771e-08 3.1736e-08
FIN 7 6.1771e-08 2.1591e-11 -- success
Gas Temperature = 1400
Gas Pressure = 101325
Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
@ -42,28 +42,28 @@ Gas Phase: gas (0)
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1400
Bulk Pressure = 101325
Bulk Temperature = 1.4e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293.065 1 1.631252e-04
Bulk Weight Growth Rate = 0.0019593 kg/m^2/s
Bulk Growth Rate = 5.56619e-07 m/s
Bulk Growth Rate = 2003.83 microns / hour
Density of bulk phase = 3520 kg / m^3
Bulk Weight Growth Rate = 0.00196 kg/m^2/s
Bulk Growth Rate = 5.57e-07 m/s
Bulk Growth Rate = 2e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1400
Surface Pressure = 101325
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0184677 -8.195044e-20
1 Csoot-H 0.981532 -8.195044e-20
Sum of coverages = 1
Surface Phase: soot_interface (11)
Surface Temperature = 1400
Surface Pressure = 101325
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0184677 -8.195044e-20
1 Csoot-H 0.981532 -8.195044e-20
@ -83,8 +83,8 @@ Sum of coverages = 1
1 5.3218e+03 2.7005e+03
2 8.9325e-06 4.0777e-06
FIN 2 8.9325e-06 6.2499e-12 -- success
Gas Temperature = 1400
Gas Pressure = 101325
Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
@ -101,28 +101,28 @@ Gas Phase: gas (0)
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1400
Bulk Pressure = 101325
Bulk Temperature = 1.4e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293.065 1 1.646335e-04
Bulk Weight Growth Rate = 0.00197741 kg/m^2/s
Bulk Growth Rate = 5.61765e-07 m/s
Bulk Growth Rate = 2022.35 microns / hour
Density of bulk phase = 3520 kg / m^3
Bulk Weight Growth Rate = 0.00198 kg/m^2/s
Bulk Growth Rate = 5.62e-07 m/s
Bulk Growth Rate = 2.02e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1400
Surface Pressure = 101325
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0186067 -3.028143e-20
1 Csoot-H 0.981393 -3.028143e-20
Sum of coverages = 1
Surface Phase: soot_interface (11)
Surface Temperature = 1400
Surface Pressure = 101325
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0186067 -3.028143e-20
1 Csoot-H 0.981393 -3.028143e-20
@ -142,8 +142,8 @@ Sum of coverages = 1
1 2.1569e+05 9.5571e+04
2 1.7622e-04 2.0108e-04
FIN 2 1.7622e-04 2.0335e-10 -- success
Gas Temperature = 1495
Gas Pressure = 101325
Gas Temperature = 1.5e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
@ -160,28 +160,28 @@ Gas Phase: gas (0)
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1495
Bulk Pressure = 101325
Bulk Temperature = 1.5e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293.065 1 2.619849e-04
Bulk Weight Growth Rate = 0.0031467 kg/m^2/s
Bulk Growth Rate = 8.93949e-07 m/s
Bulk Growth Rate = 3218.22 microns / hour
Density of bulk phase = 3520 kg / m^3
Bulk Weight Growth Rate = 0.00315 kg/m^2/s
Bulk Growth Rate = 8.94e-07 m/s
Bulk Growth Rate = 3.22e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1495
Surface Pressure = 101325
Surface Temperature = 1.5e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0242812 -8.654559e-19
1 Csoot-H 0.975719 -8.654559e-19
Sum of coverages = 1
Surface Phase: soot_interface (11)
Surface Temperature = 1495
Surface Pressure = 101325
Surface Temperature = 1.5e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0242812 -8.654559e-19
1 Csoot-H 0.975719 -8.654559e-19
@ -200,8 +200,8 @@ Sum of coverages = 1
-----------------------------------------------------------------------------------
1 2.8875e-10 1.4792e-10
FIN 1 2.8875e-10 1.2324e-10 -- success
Gas Temperature = 1495
Gas Pressure = 101325
Gas Temperature = 1.5e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (0)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
@ -218,28 +218,28 @@ Gas Phase: gas (0)
Sum of gas mole fractions= 1
Bulk Phase: soot (10)
Bulk Temperature = 1495
Bulk Pressure = 101325
Bulk Temperature = 1.5e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293.065 1 2.619849e-04
Bulk Weight Growth Rate = 0.0031467 kg/m^2/s
Bulk Growth Rate = 8.93949e-07 m/s
Bulk Growth Rate = 3218.22 microns / hour
Density of bulk phase = 3520 kg / m^3
Bulk Weight Growth Rate = 0.00315 kg/m^2/s
Bulk Growth Rate = 8.94e-07 m/s
Bulk Growth Rate = 3.22e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (11)
Surface Temperature = 1495
Surface Pressure = 101325
Surface Temperature = 1.5e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0242812 2.964615e-20
1 Csoot-H 0.975719 2.964615e-20
Sum of coverages = 1
Surface Phase: soot_interface (11)
Surface Temperature = 1495
Surface Pressure = 101325
Surface Temperature = 1.5e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0242812 2.964615e-20
1 Csoot-H 0.975719 2.964615e-20