From 53c4da6ddc452089874ef327b86b22d29d75fdac Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 16 Oct 2008 01:00:55 +0000 Subject: [PATCH] Changed the test so that it is more lenient on multiple platforms. --- test_problems/surfSolverTest/.cvsignore | 8 +- .../surfSolverTest/results2_blessed.txt | 44 ++ .../surfSolverTest/results_blessed.txt | 148 ++++ test_problems/surfSolverTest/runtest | 31 +- test_problems/surfSolverTest/runtest2 | 34 +- .../surfSolverTest/surfaceSolver.cpp | 594 ++++++++-------- .../surfSolverTest/surfaceSolver2.cpp | 658 +++++++++--------- .../surfSolverTest/surfaceSolver2_blessed.out | 96 +-- 8 files changed, 944 insertions(+), 669 deletions(-) create mode 100644 test_problems/surfSolverTest/results2_blessed.txt create mode 100644 test_problems/surfSolverTest/results_blessed.txt diff --git a/test_problems/surfSolverTest/.cvsignore b/test_problems/surfSolverTest/.cvsignore index d82c900f0..e1ed9cfbf 100644 --- a/test_problems/surfSolverTest/.cvsignore +++ b/test_problems/surfSolverTest/.cvsignore @@ -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 diff --git a/test_problems/surfSolverTest/results2_blessed.txt b/test_problems/surfSolverTest/results2_blessed.txt new file mode 100644 index 000000000..12b175a33 --- /dev/null +++ b/test_problems/surfSolverTest/results2_blessed.txt @@ -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 diff --git a/test_problems/surfSolverTest/results_blessed.txt b/test_problems/surfSolverTest/results_blessed.txt new file mode 100644 index 000000000..c53628068 --- /dev/null +++ b/test_problems/surfSolverTest/results_blessed.txt @@ -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 diff --git a/test_problems/surfSolverTest/runtest b/test_problems/surfSolverTest/runtest index 3d77a1e71..1f1ed8cfd 100755 --- a/test_problems/surfSolverTest/runtest +++ b/test_problems/surfSolverTest/runtest @@ -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 + + diff --git a/test_problems/surfSolverTest/runtest2 b/test_problems/surfSolverTest/runtest2 index cff409b21..76ca39fca 100755 --- a/test_problems/surfSolverTest/runtest2 +++ b/test_problems/surfSolverTest/runtest2 @@ -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 diff --git a/test_problems/surfSolverTest/surfaceSolver.cpp b/test_problems/surfSolverTest/surfaceSolver.cpp index a8e985ed1..3474a45e5 100644 --- a/test_problems/surfSolverTest/surfaceSolver.cpp +++ b/test_problems/surfSolverTest/surfaceSolver.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #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(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(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 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 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; } /***********************************************************/ diff --git a/test_problems/surfSolverTest/surfaceSolver2.cpp b/test_problems/surfSolverTest/surfaceSolver2.cpp index c2d711a73..4dbaff689 100644 --- a/test_problems/surfSolverTest/surfaceSolver2.cpp +++ b/test_problems/surfSolverTest/surfaceSolver2.cpp @@ -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(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(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 phaseList; - phaseList.push_back(gasTP); - phaseList.push_back(bulkPhaseTP); - phaseList.push_back(surfPhaseTP); + vector 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 phaseList2; - phaseList2.push_back(gasTP); - phaseList2.push_back(bulkPhaseTP); - phaseList2.push_back(surfPhaseTP2); + vector 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 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 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; } /***********************************************************/ diff --git a/test_problems/surfSolverTest/surfaceSolver2_blessed.out b/test_problems/surfSolverTest/surfaceSolver2_blessed.out index b7eaec50b..605b55259 100644 --- a/test_problems/surfSolverTest/surfaceSolver2_blessed.out +++ b/test_problems/surfSolverTest/surfaceSolver2_blessed.out @@ -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