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

*** empty log message ***

Browse source
This commit is contained in:
Dave Goodwin 2003-08-26 03:25:39 +00:00
parent 3b4ea774b5
commit 5d1b37c7aa
22 changed files with 316 additions and 143 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

114
Cantera/src/ChemEquil.cpp
View file

@ -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);

5
Cantera/src/ChemEquil.h
View file

@ -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;
};

22
Cantera/src/IdealGasPhase.cpp
View file

@ -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) {

2
Cantera/src/Makefile.in
View file

@ -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)

23
Cantera/src/converters/CKParser.cpp
View file

@ -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<string>& names,
bool CKParser::readThermoSection(vector<string>& 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],"<ALL>"));
if (getAllSpecies) names.clear();
map<string, int> 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;

2
Cantera/src/converters/CKParser.h
View file

@ -53,7 +53,7 @@ namespace ckr {
bool readElementSection(elementList& elements);
bool readSpeciesSection(speciesList& species);
bool readThermoSection(const vector<string>& names,
bool readThermoSection(vector<string>& names,
speciesTable& speciesData, vector_fp& temp,
int& optionFlag, ostream& log);
bool readReactionSection(const vector<string>& speciesNames,

25
Cantera/src/converters/CKReader.cpp
View file

@ -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<string> undef;
bool allsp = (speciesSymbols[0] == "<ALL>");
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 == "<empty>") {
@ -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;

33
Cantera/src/converters/ck2ct.cpp
View file

@ -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;

1
Cantera/src/ct2ctml.cpp
View file

@ -14,6 +14,7 @@
#include <string>
#include <stdlib.h>
#include "ctml.h"
#include "../../include/pypath.h"
using namespace Cantera;

21
Cantera/src/misc.cpp
View file

@ -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

20
Cantera/src/oneD/Inlet1D.h
View file

@ -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;
};
}

9
Cantera/src/oneD/MultiJac.cpp
View file

@ -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));
//}
}
}
}

9
Cantera/src/oneD/MultiNewton.cpp
View file

@ -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");
}
}

2
Cantera/src/oneD/OneDim.cpp
View file

@ -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),

10
Cantera/src/oneD/Sim1D.cpp
View file

@ -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);
}
}
}

2
Cantera/src/oneD/Sim1D.h
View file

@ -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);

2
Cantera/src/oneD/StFlow.cpp
View file

@ -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;

1
Cantera/src/oneD/StFlow.h
View file

@ -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.

75
Cantera/src/oneD/boundaries1D.cpp
View file

@ -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) {

75
Cantera/src/oneD/refine.cpp
View file

@ -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");
}
}

4
Cantera/src/oneD/refine.h
View file

@ -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);

2
Cantera/src/zeroD/Reactor.cpp
View file

@ -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();