From 5d1b37c7aa0cd5d67264510af68d73da7ee176f4 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 26 Aug 2003 03:25:39 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/src/ChemEquil.cpp | 114 ++++++++++++++++++++++------ Cantera/src/ChemEquil.h | 5 +- Cantera/src/IdealGasPhase.cpp | 22 ++---- Cantera/src/Makefile.in | 2 +- Cantera/src/converters/CKParser.cpp | 23 ++++-- Cantera/src/converters/CKParser.h | 2 +- Cantera/src/converters/CKReader.cpp | 25 +++++- Cantera/src/converters/ck2ct.cpp | 33 +++----- Cantera/src/ct2ctml.cpp | 1 + Cantera/src/misc.cpp | 21 ++++- Cantera/src/oneD/Inlet1D.h | 20 ++++- Cantera/src/oneD/MultiJac.cpp | 9 +++ Cantera/src/oneD/MultiNewton.cpp | 9 ++- Cantera/src/oneD/OneDim.cpp | 2 +- Cantera/src/oneD/Sim1D.cpp | 10 ++- Cantera/src/oneD/Sim1D.h | 2 +- Cantera/src/oneD/StFlow.cpp | 2 +- Cantera/src/oneD/StFlow.h | 1 - Cantera/src/oneD/boundaries1D.cpp | 75 +++++++++++++----- Cantera/src/oneD/refine.cpp | 75 +++++++++++------- Cantera/src/oneD/refine.h | 4 +- Cantera/src/zeroD/Reactor.cpp | 2 +- 22 files changed, 316 insertions(+), 143 deletions(-) diff --git a/Cantera/src/ChemEquil.cpp b/Cantera/src/ChemEquil.cpp index 5a4327db1..eccc574a4 100755 --- a/Cantera/src/ChemEquil.cpp +++ b/Cantera/src/ChemEquil.cpp @@ -52,7 +52,8 @@ namespace Cantera { /// Default Constructor. - ChemEquil::ChemEquil() : m_skip(-1), m_p1(0), m_p2(0), m_p0(OneAtm) + ChemEquil::ChemEquil() : m_skip(-1), m_p1(0), m_p2(0), m_p0(OneAtm), m_eloc(-1), + m_abscharge(Tiny) {} @@ -84,18 +85,56 @@ namespace Cantera { // allocate space in internal work arrays m_molefractions.resize(m_kk); - m_lambda.resize(m_mm, -10.0); + m_lambda.resize(m_mm, -100.0); m_elementmolefracs.resize(m_mm); m_comp.resize(m_mm * m_kk); m_jwork1.resize(m_mm+2); m_jwork2.resize(m_mm+2); m_startSoln.resize(m_mm+1); m_grt.resize(m_kk); + m_mu_RT.resize(m_kk); // set up elemental composition matrix - for (int k = 0; k < m_kk; k++) - for (int m = 0; m < m_mm; m++) + int m, k, mneg = -1; + doublereal na, ewt; + for (m = 0; m < m_mm; m++) { + for (k = 0; k < m_kk; k++) { + na = m_phase->nAtoms(k,m); + if (na < 0.0) { + if (mneg >= 0 && mneg != m) + throw CanteraError("ChemEquil::initialize", + "negative atom numbers allowed for only one element"); + mneg = m; + ewt = m_phase->atomicWeight(m); + if (ewt > 1.0e-3) + writelog(string("WARNING: species "+m_phase->speciesName(k) + +" has "+fp2str(m_phase->nAtoms(k,m))+" atoms of element " + +m_phase->elementName(m)+", but this element is not an electron.\n")); + } + } + } + m_eloc = mneg; +// nneg = 0.0; +// if (nneg > 0) { +// for (k = 0; k < m_kk; k++) { +// m_comp[k*m_mm + mneg] = m_phase->nAtoms(k,mneg); +// for (m = 0; m < m_mm; m++) { +// if (m != mneg) { +// m_comp[k*m_mm + m] = m_phase->nAtoms(k,m); +// m_comp[k*m_mm + mneg] += m_phase->nAtoms(k,m)*(nneg + 1); +// } +// } +// cout << m_phase->speciesName(k) << " "; +// for (m = 0; m < m_mm; m++) cout << m_comp[k*m_mm + m] << " "; +// cout << endl; +// } +// } +// else { + for (k = 0; k < m_kk; k++) { + for (m = 0; m < m_mm; m++) { m_comp[k*m_mm + m] = m_phase->nAtoms(k,m); + } + } } @@ -111,9 +150,14 @@ namespace Cantera { void ChemEquil::setToEquilState(thermo_t& s, const vector_fp& lambda_RT, doublereal t) { + fill(m_mu_RT.begin(), m_mu_RT.end(), 0.0); + for (int k = 0; k < m_kk; k++) + for (int m = 0; m < m_mm; m++) + m_mu_RT[k] += lambda_RT[m]*nAtoms(k,m); + // set the temperature - m_phase->setTemperature(t); - s.setToEquilState(lambda_RT.begin()); + s.setTemperature(t); + s.setToEquilState(m_mu_RT.begin()); update(s); } @@ -128,22 +172,24 @@ namespace Cantera { // elemental mole fractions doublereal sum = 0.0; - int m; + int m, k; for (m = 0; m < m_mm; m++) { m_elementmolefracs[m] = 0.0; - for (int k = 0; k < m_kk; k++) { + for (k = 0; k < m_kk; k++) { m_elementmolefracs[m] += nAtoms(k,m) * m_molefractions[k]; - if (nAtoms(k,m) < 0.0) { - throw CanteraError("update","negative nAtoms"); - } + //if (nAtoms(k,m) < 0.0) { + // throw CanteraError("update","negative nAtoms"); + //} if (m_molefractions[k] < 0.0) { throw CanteraError("update", "negative mole fraction for "+m_phase->speciesName(k)+ ": "+fp2str(m_molefractions[k])); } } + //cout << "update: " << m << " " << m_elementmolefracs[m] << endl; sum += m_elementmolefracs[m]; } + // normalize the element mole fractions for (m = 0; m < m_mm; m++) m_elementmolefracs[m] /= sum; } @@ -173,11 +219,11 @@ namespace Cantera { // first column contains fixed element moles for (m = 0; m < mm; m++) { aa(m+1,0) = elementMoles[m]; - if (elementMoles[m] < 0.0) { - throw CanteraError("setInitialMoles", - "negative element moles for " - +m_phase->elementName(m)+": "+fp2str(elementMoles[m])); - } + //if (elementMoles[m] < 0.0) { + // throw CanteraError("setInitialMoles", + // "negative element moles for " + // +m_phase->elementName(m)+": "+fp2str(elementMoles[m])); + // } } @@ -215,6 +261,7 @@ namespace Cantera { for (int k = 0; k < kksp; k++) { if (ip == ksp) { m_molefractions[k] = aa(n+1, 0); + //cout << "initial " << m_phase->speciesName(k) << " " << m_molefractions[k] << endl; } ksp++; } @@ -351,15 +398,15 @@ namespace Cantera { int nvar = mm + 1; DenseMatrix jac(nvar, nvar); // jacobian - vector_fp x(nvar, -10.0); // solution vector + vector_fp x(nvar, -100.0); // solution vector vector_fp res_trial(nvar); vector_fp elementMol(mm, 0.0); double perturb; for (m = 0; m < mm; m++) { + if (m_skip < 0 && elMoles[m] > 0.0 ) m_skip = m; perturb = Cutoff*(1.0 + rand()); elementMol[m] = elMoles[m] + perturb; - if (m_skip < 0 && elMoles[m] > 0.0 ) m_skip = m; } update(s); @@ -427,11 +474,11 @@ namespace Cantera { setInitialMoles(s, elementMol); - for (int ii = 0; ii < m_mm; ii++) x[ii] = -10.0; - try { + for (int ii = 0; ii < m_mm; ii++) x[ii] = -100.0; + //try { estimateElementPotentials(s, x); - } - catch (CanteraError) { ; } + //} + //catch (CanteraError) { ; } x[m_mm] = log(m_phase->temperature()); @@ -441,7 +488,8 @@ namespace Cantera { for (m = 0; m < mm; m++) { above[m] = 200.0; below[m] = -2000.0; - if (elMoles[m] < Cutoff) x[m] = -1000.0; + if (elMoles[m] < Cutoff && m != m_eloc) x[m] = -1000.0; + //if (m == m_eloc) x[m] = -10.0; } above[mm] = log(1.e4); below[mm] = log(10.0); @@ -458,6 +506,12 @@ namespace Cantera { next: + if (m_eloc >= 0) { + m_abscharge = 0.0; + int k; + for (k = 0; k < m_kk; k++) m_abscharge += fabs(m_phase->charge(k)*m_molefractions[k]); + } + iter++; equilResidual(s, x, elMoles, res_trial, XY, xval, yval); f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin()); @@ -515,6 +569,7 @@ namespace Cantera { yy = m_p2->value(s); deltax = (xx - xval)/xval; deltay = (yy - yval)/yval; + if (absmax(res_trial.begin(), res_trial.end()) < options.relTolerance && fabs(deltax) < options.relTolerance && fabs(deltay) < options.relTolerance) { @@ -639,11 +694,20 @@ namespace Cantera { for (n=0; n < m_mm; n++) { // drive element potential for absent elements to -1000 - if (elmtotal[n] < Cutoff) + if (elmtotal[n] < Cutoff && n != m_eloc) resid[n] = x[n] + 1000.0; else resid[n] = log( (1.0 + elmtotal[n]) / (1.0 + elm[n]) ); } + if (m_eloc >= 0) { + doublereal chrg, sumnet = 0.0, sumabs = 0.0; + for (int k = 0; k < m_kk; k++) { + chrg = m_molefractions[k]*m_phase->charge(k); + sumnet += chrg; + sumabs += fabs(chrg); + } + resid[m_eloc] = sumnet/m_abscharge; // log((1.0 + sumnet/sumabs)); + } xx = m_p1->value(mix); yy = m_p2->value(mix); resid[m_mm] = xx/xval - 1.0; @@ -665,7 +729,7 @@ namespace Cantera { int n, m; doublereal rdx, dx, xsave; - doublereal atol = 1.e-7; + doublereal atol = 1.e-10; equilResidual(mix, x, elmols, r0, XY, xval, yval); diff --git a/Cantera/src/ChemEquil.h b/Cantera/src/ChemEquil.h index 3b02f6016..2d5c5c272 100755 --- a/Cantera/src/ChemEquil.h +++ b/Cantera/src/ChemEquil.h @@ -39,7 +39,7 @@ namespace Cantera { */ class EquilOpt { public: - EquilOpt() : relTolerance(1.e-9), maxIterations(1000), iterations(0), + EquilOpt() : relTolerance(1.e-10), maxIterations(1000), iterations(0), maxStepSize(10.0), propertyPair(TP), contin(false) {} doublereal relTolerance; ///< Relative tolerance @@ -136,11 +136,14 @@ namespace Cantera { vector_fp m_comp; doublereal m_temp, m_dens; doublereal m_p0; + int m_eloc; + doublereal m_abscharge; doublereal m_startTemp, m_startDens; vector_fp m_startSoln; vector_fp m_grt; + vector_fp m_mu_RT; }; diff --git a/Cantera/src/IdealGasPhase.cpp b/Cantera/src/IdealGasPhase.cpp index 7b4b5c57d..430a895af 100644 --- a/Cantera/src/IdealGasPhase.cpp +++ b/Cantera/src/IdealGasPhase.cpp @@ -51,35 +51,25 @@ namespace Cantera { /** - * Set mixture to an equilibrium state consistent with specified - * element potentials and temperature. - * - * @param lambda_RT vector of non-dimensional element potentials - * \f[ \lambda_m/RT \f]. - * @param t temperature in K. - * @param work. Temporary work space. Must be dimensioned at least - * as large as the number of species. - * + * Set mixture to an equilibrium state consistent with specified + * chemical potentials and temperature. */ - void IdealGasPhase::setToEquilState(const doublereal* lambda_RT) + void IdealGasPhase::setToEquilState(const doublereal* mu_RT) { const array_fp& grt = gibbs_RT(); - // set the pressure and composition to be consistent with - // the temperature, doublereal pres = 0.0; for (int k = 0; k < m_kk; k++) { - m_pp[k] = -grt[k]; - for (int m = 0; m < m_mm; m++) { - m_pp[k] += nAtoms(k,m)*lambda_RT[m]; - } + m_pp[k] = -grt[k] + mu_RT[k]; m_pp[k] = m_p0 * exp(m_pp[k]); pres += m_pp[k]; + //cout <<"setToEquil: " << k << " " << grt[k] << " " << mu_RT[k] << endl; } // set state setState_PX(pres, m_pp.begin()); } + void IdealGasPhase::_updateThermo() const { doublereal tnow = temperature(); if (m_tlast != tnow) { diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in index 9ef79f5fe..52d8b004e 100755 --- a/Cantera/src/Makefile.in +++ b/Cantera/src/Makefile.in @@ -51,7 +51,7 @@ SOLVERS = CVode.o BandMatrix.o FLOW1D = $(KINETICS) $(SOLVERS) -EVERYTHING = $(KINETICS) $(HETEROKIN) $(CK) $(TRANSPORT) $(REACTOR) $(RPATH) \ +EVERYTHING = $(KINETICS) $(HETEROKIN) $(EQUIL) $(CK) $(TRANSPORT) $(REACTOR) $(RPATH) \ $(SOLVERS) $(FLOW1D) diff --git a/Cantera/src/converters/CKParser.cpp b/Cantera/src/converters/CKParser.cpp index bde7defd5..99bf99334 100755 --- a/Cantera/src/converters/CKParser.cpp +++ b/Cantera/src/converters/CKParser.cpp @@ -6,7 +6,10 @@ // Copyright 2001 California Institute of Technology // // $Log$ -// Revision 1.2 2003-07-30 20:56:57 dggoodwin +// Revision 1.3 2003-08-26 03:25:41 dggoodwin +// *** empty log message *** +// +// Revision 1.2 2003/07/30 20:56:57 dggoodwin // *** empty log message *** // // Revision 1.1.1.1 2003/04/14 17:57:52 dggoodwin @@ -161,10 +164,10 @@ namespace ckr { throw CK_SyntaxError(log, "error reading Tmin, Tmid, or Tmax"); } - if (tmin > tmid || tmid > tmax) { - throw CK_SyntaxError(log, - "condition Tmin <= Tmid <= Tmax violated"); - } + //if (tmin > tmid || tmid > tmax) { + // throw CK_SyntaxError(log, + // "condition Tmin <= Tmid <= Tmax violated"); + //} } @@ -526,7 +529,7 @@ next: * */ - bool CKParser::readThermoSection(const vector& names, + bool CKParser::readThermoSection(vector& names, speciesTable& species, vector_fp& temp, int& optionFlag, ostream& log) { string s; @@ -603,6 +606,7 @@ next: // now read in all species records that have names in list 'names' bool getAllSpecies = (nsp > 0 && match(names[0],"")); + if (getAllSpecies) names.clear(); map dup; // used to check for duplicate THERMO records bool already_read; @@ -651,7 +655,12 @@ next: writeSpeciesData(log, spec); } checkTemps(log, spec.tlow, spec.tmid, spec.thigh); - nsp--; + if (getAllSpecies) { + names.push_back(spec.name); + nsp = names.size(); + } + else + nsp--; } } return true; diff --git a/Cantera/src/converters/CKParser.h b/Cantera/src/converters/CKParser.h index 4705c8f00..dcc92eca1 100755 --- a/Cantera/src/converters/CKParser.h +++ b/Cantera/src/converters/CKParser.h @@ -53,7 +53,7 @@ namespace ckr { bool readElementSection(elementList& elements); bool readSpeciesSection(speciesList& species); - bool readThermoSection(const vector& names, + bool readThermoSection(vector& names, speciesTable& speciesData, vector_fp& temp, int& optionFlag, ostream& log); bool readReactionSection(const vector& speciesNames, diff --git a/Cantera/src/converters/CKReader.cpp b/Cantera/src/converters/CKReader.cpp index 8c5fda25d..a5bdb0860 100755 --- a/Cantera/src/converters/CKReader.cpp +++ b/Cantera/src/converters/CKReader.cpp @@ -160,16 +160,23 @@ bool CKReader::read(const string& inputFile, const string& thermoDatabase, bool hasthermo = parser.advanceToKeyword("THERM","REAC"); - int optionFlag = 0; + int k, optionFlag = 0; int undefined = species.size(); string nm; vector undef; - + bool allsp = (speciesSymbols[0] == ""); if (hasthermo && parser.readThermoSection(speciesSymbols, speciesData, temp, optionFlag, log)) { - + if (allsp) { + nsp = speciesData.size(); + for (k = 0; k < nsp; k++) { + Species s; + s.name = speciesSymbols[k]; + species.push_back(s); + } + } undefined = 0; - for (int k = 0; k < nsp; k++) { + for (k = 0; k < nsp; k++) { nm = species[k].name; species[k] = speciesData[species[k].name]; if (species[k].name == "") { @@ -204,6 +211,16 @@ bool CKReader::read(const string& inputFile, const string& thermoDatabase, int dbflag = HasTempRange; thermoReader.readThermoSection(undef, speciesData, temp, dbflag, log); undefined = 0; + if (allsp) { + species.clear(); + nsp = speciesData.size(); + for (k = 0; k < nsp; k++) { + Species s; + s.name = undef[k]; + cout << k << " " << s.name << endl; + species.push_back(s); + } + } for (int k = 0; k < nsp; k++) { if (species[k].valid == 0) { nm = species[k].name; diff --git a/Cantera/src/converters/ck2ct.cpp b/Cantera/src/converters/ck2ct.cpp index 2b4cf6e14..1fa0f98a5 100644 --- a/Cantera/src/converters/ck2ct.cpp +++ b/Cantera/src/converters/ck2ct.cpp @@ -217,15 +217,13 @@ namespace pip { if (nm == "E") charge = -sp.elements[m].number; } - /* - * Add the child element, atomArray, to the species xml node. - */ printf(" atoms = \"%s\",\n", str.c_str()); addNASA(sp.lowCoeffs, sp.highCoeffs, sp.tlow, sp.tmid, sp.thigh); if (_with_transport) addTransportParams(sp.name); + if (sp.id != "") printf(",\n note = \"%s\"", sp.id.c_str()); printf("\n )\n"); } @@ -305,19 +303,11 @@ namespace pip { * mechanism. The node will be filled up with the description * of the mechanism. This is the output to the routine. */ - void ck2ct(string idtag, ckr::CKReader& r) { + void ck2ct(string idtag, ckr::CKReader& r, bool hastransport) { popError(); doublereal version = 1.0; - //cout << "dataset(\"" << idtag << "\")" << endl; - - // cout << "\n\n"; - //writeline(); - //cout << "#\n# The default units. These will be used for dimensional quantities" << endl - // << "# with unspecified units." << endl; - //writeline(); - cout << "units(length = \"cm\", time = \"s\", quantity = \"mol\", "; string e_unit; int eunit = r.units.ActEnergy; @@ -335,13 +325,6 @@ namespace pip { e_unit = "eV"; cout << "act_energy = " << "\"" << e_unit << "\")\n\n"; - //cout << "\n\n\n"; - //writeline(); - //cout << "#\n# The phase definition. This specifies an ideal gas mixture that" << endl - // << "# includes all species and reactions defined in this file." - // << "\n#\n"; - //writeline(); - printf("\nideal_gas(name = \"%s\",\n",idtag.c_str()); @@ -367,7 +350,12 @@ namespace pip { if ((i+1) % 10 == 0) spnames += "\n "; } printf(" species = \"\"\"%s\"\"\",\n", spnames.c_str()); - printf(" reactions = \"all\""); + printf(" reactions = \"all\",\n"); + if (hastransport) { + printf(" transport = \"Mix\",\n"); + } + printf(" initial_state = state(temperature = 300.0,\n"); + printf(" pressure = OneAtm)"); cout << " )" << endl; cout << "\n\n\n"; @@ -423,7 +411,7 @@ namespace pip { const char* tr_file, const char* id_tag) { ckr::CKReader r; - r.validate = false; + r.validate = true; //int i=1; string infile = string(in_file); @@ -456,7 +444,8 @@ namespace pip { cout << "# Transport data from file "+trfile+".\n" << endl; getTransportData(trfile); } - ck2ct(idtag, r); + bool hastransport = (trfile != ""); + ck2ct(idtag, r, hastransport); } catch (CanteraError) { return -1; diff --git a/Cantera/src/ct2ctml.cpp b/Cantera/src/ct2ctml.cpp index 0d99f23d9..958afb327 100644 --- a/Cantera/src/ct2ctml.cpp +++ b/Cantera/src/ct2ctml.cpp @@ -14,6 +14,7 @@ #include #include #include "ctml.h" +#include "../../include/pypath.h" using namespace Cantera; diff --git a/Cantera/src/misc.cpp b/Cantera/src/misc.cpp index f4ef9d2d5..6133500ac 100755 --- a/Cantera/src/misc.cpp +++ b/Cantera/src/misc.cpp @@ -93,14 +93,26 @@ namespace Cantera { } if (ext != ".xml" && ext != ".ctml") { ctml::ct2ctml(path.c_str()); - ff = path.substr(0,idot) + ".xml"; + string::size_type islash = path.rfind('/'); + if (islash != string::npos) + ff = string("./")+path.substr(islash+1,idot-islash - 1) + ".xml"; + else + ff = string("./")+path.substr(0,idot) + ".xml"; + writelog("ff = "+ff+"\n"); } else { ff = path; } ifstream s(ff.c_str()); XML_Node* x = new XML_Node("doc"); - x->build(s); + if (s) { + x->build(s); + __app->xmlfiles[file] = x; + __app->xmlfiles[ff] = x; + } + else + throw CanteraError("get_XML_File","cannot open "+ff+" for reading."); + /* * Add the built XML Tree to the map, xmlfiles. * It stores the pointer to the tree, with the @@ -109,7 +121,6 @@ namespace Cantera { * HKM Note: shouldn't the key be the full pathname of * the file, i.e., ff? */ - __app->xmlfiles[file] = x; } return __app->xmlfiles[file]; } @@ -290,11 +301,15 @@ namespace Cantera { inname = ""; for (i = 0; i < nd; i++) { inname = dirs[i] + "/" + name; + writelog("looking for file "+inname+"\n"); ifstream fin(inname.c_str()); if (fin) { + writelog("found it\n"); fin.close(); return inname; } + else + writelog("not found.\n"); } string msg; msg = "\nInput file " + name diff --git a/Cantera/src/oneD/Inlet1D.h b/Cantera/src/oneD/Inlet1D.h index e0f9e1a53..2ddc70778 100644 --- a/Cantera/src/oneD/Inlet1D.h +++ b/Cantera/src/oneD/Inlet1D.h @@ -72,7 +72,7 @@ namespace Cantera { virtual doublereal mdot() {return m_mdot;} virtual void _getInitialSoln(doublereal* x) { - cout << "Bdry1D::_getInitialSoln called! " << m_index << endl; + writelog("Bdry1D::_getInitialSoln called!\n"); } protected: @@ -112,6 +112,7 @@ namespace Cantera { Inlet1D() : Bdry1D(), m_V0(0.0), m_nsp(0), m_flow(0) { m_type = cInletType; m_xstr = ""; + writelog("Inlet1D constructor\n"); } virtual ~Inlet1D(){} @@ -217,6 +218,10 @@ namespace Cantera { virtual void _finalize(const doublereal* x) { ; //m_temp = x[0]; } + virtual void _getInitialSoln(doublereal* x) { + x[0] = m_temp; + } + protected: }; @@ -245,6 +250,9 @@ namespace Cantera { virtual void _finalize(const doublereal* x) { ; //m_temp = x[0]; } + virtual void _getInitialSoln(doublereal* x) { + x[0] = m_temp; + } protected: }; @@ -262,6 +270,7 @@ namespace Cantera { Surf1D() : Bdry1D() { m_type = cSurfType; + writelog("Surf1D constructor\n"); } virtual ~Surf1D(){} @@ -314,12 +323,15 @@ namespace Cantera { } void setKineticsMgr(InterfaceKinetics* kin) { + m_kin = kin; m_surfindex = kin->surfacePhaseIndex(); m_sphase = (SurfPhase*)&kin->thermo(m_surfindex); m_nsp = m_sphase->nSpecies(); m_enabled = true; } + void enableCoverageEquations(bool docov) { m_enabled = docov; } + virtual ~ReactingSurf1D(){} virtual string componentName(int n) const; @@ -334,11 +346,12 @@ namespace Cantera { virtual void _getInitialSoln(doublereal* x) { x[0] = m_temp; + //m_kin->advanceCoverages(1.0); m_sphase->getCoverages(x+1); } virtual void _finalize(const doublereal* x) { - ; //m_temp = x[0]; + copy(x+1,x+1+m_nsp,m_fixed_cov.begin()); } virtual void showSolution(const doublereal* x) { @@ -360,8 +373,9 @@ namespace Cantera { SurfPhase* m_sphase; int m_surfindex, m_nsp; bool m_enabled; - vector_fp m_fixed_cov; vector_fp m_work; + vector_fp m_fixed_cov; + int dum; }; } diff --git a/Cantera/src/oneD/MultiJac.cpp b/Cantera/src/oneD/MultiJac.cpp index 2687b3838..ff8b3ad91 100644 --- a/Cantera/src/oneD/MultiJac.cpp +++ b/Cantera/src/oneD/MultiJac.cpp @@ -24,6 +24,11 @@ using namespace std; #include "MultiJac.h" +static int isnan(double x) { + if (x > 1.0e300) return 1; + return 0; +} + namespace Cantera { MultiJac::MultiJac(OneDim& r) @@ -92,6 +97,10 @@ namespace Cantera { for (m = 0; m < mv; m++) { value(m+iloc,ipt) = (m_r1[m+iloc] - resid0[m+iloc])*rdx; + //if (isnan(value(m+iloc,ipt))) { + // throw CanteraError("eval","isnan: "+fp2str(m_r1[m+iloc]) + // +" "+fp2str(resid0[m+iloc])+" "+fp2str(rdx)); + //} } } } diff --git a/Cantera/src/oneD/MultiNewton.cpp b/Cantera/src/oneD/MultiNewton.cpp index 762a3a3da..515930d80 100644 --- a/Cantera/src/oneD/MultiNewton.cpp +++ b/Cantera/src/oneD/MultiNewton.cpp @@ -120,8 +120,10 @@ namespace Cantera { for (n = 0; n < sz; n++) { step[n] = -step[n]; } - jac.solve(sz, step, step); - + try { + jac.solve(sz, step, step); + } + catch (CanteraError) { #undef DEBUG_STEP #ifdef DEBUG_STEP bool ok = false; @@ -138,6 +140,9 @@ namespace Cantera { //if (!ok) throw "not ok"; } #endif + writelog("size = "+int2str(sz)+"\n"); + throw CanteraError("step","step error"); + } } diff --git a/Cantera/src/oneD/OneDim.cpp b/Cantera/src/oneD/OneDim.cpp index 410886e54..34b122b0e 100644 --- a/Cantera/src/oneD/OneDim.cpp +++ b/Cantera/src/oneD/OneDim.cpp @@ -17,7 +17,7 @@ namespace Cantera { * Default constructor. Create an empty object. */ OneDim::OneDim() - : m_tmin(1.0e-16), m_tmax(0.1), m_tfactor(0.5), + : m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5), m_jac(0), m_newt(0), m_rdt(0.0), m_jac_ok(false), m_nd(0), m_bw(0), m_size(0), diff --git a/Cantera/src/oneD/Sim1D.cpp b/Cantera/src/oneD/Sim1D.cpp index af75eb06b..57b060495 100644 --- a/Cantera/src/oneD/Sim1D.cpp +++ b/Cantera/src/oneD/Sim1D.cpp @@ -272,9 +272,11 @@ namespace Cantera { if (loglevel > 2) showSolution(); if (refine_grid) { + //writelog("calling refine.\n"); new_points = refine(loglevel); } else { + writelog("grid refinement disabled.\n"); new_points = 0; } } @@ -327,6 +329,7 @@ namespace Cantera { zmid = 0.5*(d.grid(m) + d.grid(m+1)); znew.push_back(zmid); np++; + //writelog(string("refine: adding point at ")+fp2str(zmid)+"\n"); // for each component, linearly interpolate // the solution to this point @@ -337,6 +340,7 @@ namespace Cantera { } } else { + writelog(string("refine: discarding point at ")+fp2str(d.grid(m))+"\n"); ; // throw CanteraError("refine","keepPoint is false at m = "+int2str(m)); } } @@ -378,15 +382,15 @@ namespace Cantera { * are applied to each domain. @see Refiner::setCriteria. */ void Sim1D::setRefineCriteria(int dom, doublereal ratio, - doublereal slope, doublereal curve) { + doublereal slope, doublereal curve, doublereal prune) { if (dom >= 0) { Refiner& r = domain(dom).refiner(); - r.setCriteria(ratio, slope, curve); + r.setCriteria(ratio, slope, curve, prune); } else { for (int n = 0; n < m_nd; n++) { Refiner& r = domain(n).refiner(); - r.setCriteria(ratio, slope, curve); + r.setCriteria(ratio, slope, curve, prune); } } } diff --git a/Cantera/src/oneD/Sim1D.h b/Cantera/src/oneD/Sim1D.h index a1981cc2e..ec76f4cb5 100644 --- a/Cantera/src/oneD/Sim1D.h +++ b/Cantera/src/oneD/Sim1D.h @@ -87,7 +87,7 @@ namespace Cantera { /// Set the criteria for grid refinement. void setRefineCriteria(int dom = -1, doublereal ratio = 10.0, - doublereal slope = 0.8, doublereal curve = 0.8); + doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1); void restore(string fname, string id); diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp index 2851a8521..83dc61abb 100644 --- a/Cantera/src/oneD/StFlow.cpp +++ b/Cantera/src/oneD/StFlow.cpp @@ -431,7 +431,7 @@ namespace Cantera { doublereal sum = 0.0; for (k = 0; k < m_nsp; k++) { sum += Y(x,k,j); - rsd[index(k+4,j)] = m_flux(k,j-1); + rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); } rsd[index(4,j)] = 1.0 - sum; diag[index(4,j)] = 0; diff --git a/Cantera/src/oneD/StFlow.h b/Cantera/src/oneD/StFlow.h index 6ea148232..86461ac77 100644 --- a/Cantera/src/oneD/StFlow.h +++ b/Cantera/src/oneD/StFlow.h @@ -275,7 +275,6 @@ namespace Cantera { void setGas(const doublereal* x,int j); void setGasAtMidpoint(const doublereal* x,int j); - protected: // used to write mass fractions to plot files. diff --git a/Cantera/src/oneD/boundaries1D.cpp b/Cantera/src/oneD/boundaries1D.cpp index 917f61583..b25769ef7 100644 --- a/Cantera/src/oneD/boundaries1D.cpp +++ b/Cantera/src/oneD/boundaries1D.cpp @@ -103,7 +103,7 @@ namespace Cantera { needJacUpdate(); } } - + string Inlet1D:: componentName(int n) const { switch (n) { @@ -203,7 +203,8 @@ namespace Cantera { rb[0] += x[0]; // u for (k = 1; k < m_nsp; k++) { if (m_flow->doSpecies(k)) { - rb[4+k] += x[0]*(-xb[4+k] + m_yin[k]); + // rb[4+k] += x[0]*(-xb[4+k] + m_yin[k]); + rb[4+k] += x[0]*(m_yin[k]); } } } @@ -352,7 +353,7 @@ namespace Cantera { r[0] = x[0]; diag[0] = 0; - int nc; + int nc, k; if (m_flow_right) { nc = m_flow_right->nComponents(); @@ -361,6 +362,11 @@ namespace Cantera { db = diag + 1; rb[0] = xb[3]; rb[2] = xb[2] - xb[2 + nc]; + for (k = 4; k < nc; k++) { + //if (m_flow_right->doSpecies(k-4)) { + rb[k] = xb[k] - xb[k + nc]; + //} + } } if (m_flow_left) { @@ -370,6 +376,11 @@ namespace Cantera { db = diag - nc; rb[0] = xb[3]; rb[2] = xb[2] - xb[2 - nc]; + for (k = 5; k < nc; k++) { + // if (m_flow_left->doSpecies(k-4)) { + rb[k] = xb[k] - xb[k - nc]; + //} + } } } @@ -490,22 +501,24 @@ namespace Cantera { } void ReactingSurf1D:: - init() { + init() { + m_nv = m_nsp + 1; _init(m_nsp+1); m_fixed_cov.resize(m_nsp, 0.0); m_fixed_cov[0] = 1.0; - m_work.resize(m_kin->nTotalSpecies()); + int nt = m_kin->nTotalSpecies(); + m_work.resize(nt, 0.0); + // set bounds vector_fp lower(m_nv), upper(m_nv); lower[0] = 200.0; upper[0] = 1.e5; int n; for (n = 0; n < m_nsp; n++) { - lower[n] = -1.0e-5; - upper[n] = 2.0; + lower[n+1] = -1.0e-5; + upper[n+1] = 2.0; } setBounds(m_nv, lower.begin(), m_nv, upper.begin()); - vector_fp rtol(m_nv), atol(m_nv); for (n = 0; n < m_nv; n++) { rtol[n] = 1.0e-5; @@ -534,10 +547,13 @@ namespace Cantera { doublereal sum = 0.0; int k; for (k = 0; k < m_nsp; k++) { - m_work[k] = x[k]; - sum += x[k]; + m_work[k] = x[k+1]; + sum += x[k+1]; } + m_sphase->setTemperature(x[0]); m_sphase->setCoverages(m_work.begin()); + //m_kin->advanceCoverages(1.0); + //m_sphase->getCoverages(m_fixed_cov.begin()); // set the left gas state to the adjacent point @@ -567,22 +583,41 @@ namespace Cantera { doublereal maxx = -1.0; int imx = -1; for (k = 0; k < m_nsp; k++) { - r[k] = m_work[k + ioffset] * m_sphase->size(k) * rs0; - r[k] -= rdt*(x[k] - prevSoln(k,0)); - diag[k] = 1; - if (x[k] > maxx) { - maxx = x[k]; - imx = k; + r[k+1] = m_work[k + ioffset] * m_sphase->size(k) * rs0; + r[k+1] -= rdt*(x[k+1] - prevSoln(k+1,0)); + diag[k+1] = 1; + if (x[k+1] > maxx) { + maxx = x[k+1]; + imx = k+1; } } - r[imx] = 1.0 - sum; - diag[imx] = 0; + r[1] = 1.0 - sum; + diag[1] = 0; } else { - r[k] = x[k] - m_fixed_cov[k]; - diag[k] = 0; + for (k = 0; k < m_nsp; k++) { + r[k+1] = x[k+1] - m_fixed_cov[k]; + diag[k+1] = 0; + } } + if (m_flow_right) { + rb = r + 1; + xb = x + 1; + rb[2] = xb[2] - x[0]; // specified T + } + int nc; + if (m_flow_left) { + nc = m_flow_left->nComponents(); + const doublereal* mwleft = m_phase_left->molecularWeights().begin(); + rb =r - nc; + xb = x - nc; + rb[2] = xb[2] - x[0]; // specified T + for (int nl = 1; nl < m_left_nsp; nl++) { + rb[4+nl] += m_work[nl]*mwleft[nl]; + } + } + // gas-phase residuals // doublereal rho; // if (m_flow_left) { diff --git a/Cantera/src/oneD/refine.cpp b/Cantera/src/oneD/refine.cpp index 4779bf56f..08a1a7e68 100644 --- a/Cantera/src/oneD/refine.cpp +++ b/Cantera/src/oneD/refine.cpp @@ -33,8 +33,8 @@ namespace Cantera { Refiner::Refiner(Domain1D& domain) : - m_ratio(10.0), m_slope(0.8), m_curve(0.8), m_prune(0.1), - m_min_range(0.001), m_domain(&domain), m_npmax(300) + m_ratio(10.0), m_slope(0.8), m_curve(0.8), m_prune(-0.001), + m_min_range(0.01), m_domain(&domain), m_npmax(300) { m_nv = m_domain->nComponents(); m_active.resize(m_nv, true); @@ -45,16 +45,22 @@ namespace Cantera { int Refiner::analyze(int n, const doublereal* z, const doublereal* x) { + if (n >= m_npmax) { + writelog("max number of grid points reached ("+int2str(m_npmax)+".\n"); + return -2; + } + + if (m_domain->nPoints() <= 1) { + //writelog("can't refine a domain with 1 point: "+m_domain->id()+"\n"); + return 0; + } + m_loc.clear(); m_c.clear(); m_keep.clear(); m_keep[0] = 1; m_keep[n-1] = 1; - //m_did_analysis = false; - - - if (m_domain->nPoints() <= 1) return 0; m_nv = m_domain->nComponents(); @@ -63,10 +69,6 @@ namespace Cantera { throw CanteraError("analyze","inconsistent"); - if (n >= m_npmax) { - return -2; // throw CanteraError("analyze","max points"); - } - /** * find locations where cell size ratio is too large. */ @@ -77,24 +79,10 @@ namespace Cantera { doublereal dmax, r; vector_fp v(n), s(n-1); - dz[0] = z[1] - z[0]; - // for (j = 1; j < n-1; j++) { -// dz[j] = z[j+1] - z[j]; -// if (dz[j] > m_ratio*dz[j-1]) { -// m_loc[j] = 1; -// m_c["point "+int2str(j)] = 1; -// } -// if (dz[j] < dz[j-1]/m_ratio) { -// m_loc[j-1] = 1; -// m_c["point "+int2str(j-1)] = 1; -// } -// if (m_loc.size() + n > m_npmax) goto done; -// } - for (int i = 0; i < m_nv; i++) { if (m_active[i]) { name = m_domain->componentName(i); - + //writelog("refine: examining "+name+"\n"); // get component i at all points for (j = 0; j < n; j++) v[j] = value(x, i, j); @@ -137,8 +125,12 @@ namespace Cantera { m_keep[j+1] = 1; } else { - if (m_keep[j] == 0) m_keep[j] = -1; - if (m_keep[j+1] == 0) m_keep[j+1] = -1; + //writelog(string("r = ")+fp2str(r)+"\n"); + if (m_keep[j] == 0) { + //if (m_keep[j-1] > -1 && m_keep[j+1] > -1) + m_keep[j] = -1; + } + //if (m_keep[j+1] == 0) m_keep[j+1] = -1; } } } @@ -166,13 +158,34 @@ namespace Cantera { m_keep[j+1] = 1; } else { - if (m_keep[j+1] == 0) m_keep[j+1] = -1; + //writelog(string("r slope = ")+fp2str(r)+"\n"); + if (m_keep[j+1] == 0) { + //if (m_keep[j] > -1 && m_keep[j+2] > -1) + m_keep[j+1] = -1; + } } } } } } + + dz[0] = z[1] - z[0]; + for (j = 1; j < n-1; j++) { + dz[j] = z[j+1] - z[j]; + if (dz[j] > m_ratio*dz[j-1]) { + m_loc[j] = 1; + m_c["point "+int2str(j)] = 1; + } + if (dz[j] < dz[j-1]/m_ratio) { + m_loc[j-1] = 1; + m_c["point "+int2str(j-1)] = 1; + } + //if (m_loc.size() + n > m_npmax) goto done; + } + + + done: //m_did_analysis = true; return m_loc.size(); @@ -201,6 +214,12 @@ done: } writelog("\n"); } + else if (m_domain->nPoints() > 1) { + writelog("no new points needed in "+m_domain->id()+"\n"); + //writelog("curve = "+fp2str(m_curve)+"\n"); + //writelog("slope = "+fp2str(m_slope)+"\n"); + //writelog("prune = "+fp2str(m_prune)+"\n"); + } } diff --git a/Cantera/src/oneD/refine.h b/Cantera/src/oneD/refine.h index a9fb61a7e..97f909f54 100644 --- a/Cantera/src/oneD/refine.h +++ b/Cantera/src/oneD/refine.h @@ -15,7 +15,7 @@ namespace Cantera { void setCriteria(doublereal ratio = 10.0, doublereal slope = 0.8, doublereal curve = 0.8, - doublereal prune = 0.1) { + doublereal prune = -0.1) { m_ratio = ratio; m_slope = slope; m_curve = curve; m_prune = prune; } @@ -30,7 +30,7 @@ namespace Cantera { return m_loc.find(j) != m_loc.end(); } bool keepPoint(int j) { - return m_keep[j] != -1; // m_keep.find(j) != m_keep.end(); + return (m_keep[j] != -1); // m_keep.find(j) != m_keep.end(); } double value(const double* x, int i, int j); diff --git a/Cantera/src/zeroD/Reactor.cpp b/Cantera/src/zeroD/Reactor.cpp index f33af3194..0e888530f 100644 --- a/Cantera/src/zeroD/Reactor.cpp +++ b/Cantera/src/zeroD/Reactor.cpp @@ -112,7 +112,7 @@ namespace Cantera { m_pressure = m_thermo->pressure(); m_intEnergy = m_thermo->intEnergy_mass(); - int nt, maxnt = 0; + int nt = 0, maxnt = 0; for (int m = 0; m < m_nwalls; m++) { if (m_wall[m]->kinetics(m_lr[m])) { nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();