GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

*** empty log message ***

Browse Source
This commit is contained in:
Dave Goodwin
2003-08-13 13:56:34 +00:00
parent 4391c86e45
commit 81189197f1
17 changed files with 543 additions and 188 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
Cantera/clib/src/ctbdry.cpp
View File

@@ -12,7 +12,7 @@
// Cantera includes
#include "oneD/OneDim.h"
#include "oneD/Inlet1D.h"
#include "InterfaceKinetics.h"
#include "Cabinet.h"
#include "Storage.h"
@@ -34,6 +34,10 @@ inline ThermoPhase* _thermo(int n) {
return Storage::__storage->__thtable[n];
}
inline Kinetics* _kin(int n) {
return Storage::__storage->__ktable[n];
}
extern "C" {
int DLL_EXPORT bndry_new(int itype) {
@@ -45,6 +49,8 @@ extern "C" {
s = new Symm1D(); break;
case 3:
s = new Surf1D(); break;
case 4:
s = new ReactingSurf1D(); break;
default:
return -2;
}
@@ -114,4 +120,14 @@ extern "C" {
catch (CanteraError) {return -1;}
return 0;
}
int DLL_EXPORT surf_setkinetics(int i, int j) {
try {
ReactingSurf1D* srf = (ReactingSurf1D*)_bndry(i);
InterfaceKinetics* k = (InterfaceKinetics*)_kin(j);
srf->setKineticsMgr(k);
}
catch (CanteraError) {return -1;}
return 0;
}
}

1
Cantera/clib/src/ctbdry.h
View File

@@ -15,5 +15,6 @@ extern "C" {
double DLL_IMPORT bndry_mdot(int i);
int DLL_IMPORT bndry_setxin(int i, double* xin);
int DLL_IMPORT bndry_setxinbyname(int i, char* xin);
int DLL_IMPORT bndry_setkinetics(int i, int j);
}
#endif

17
Cantera/clib/src/ctonedim.cpp
View File

@@ -179,6 +179,14 @@ extern "C" {
catch (CanteraError) { return -1; }
}
int DLL_EXPORT reactingsurf_new() {
try {
ReactingSurf1D* i = new ReactingSurf1D();
return Cabinet<Domain1D>::cabinet()->add(i);
}
catch (CanteraError) { return -1; }
}
int DLL_EXPORT symm_new() {
try {
Symm1D* i = new Symm1D();
@@ -240,6 +248,15 @@ extern "C" {
catch (CanteraError) { return -1; }
}
int DLL_EXPORT reactingsurf_setkineticsmgr(int i, int j) {
try {
ReactingSurf1D* srf = (ReactingSurf1D*)_bdry(i);
InterfaceKinetics* k = (InterfaceKinetics*)_kinetics(j);
srf->setKineticsMgr(k);
return 0;
}
catch (CanteraError) { return -1; }
}
//------------------ stagnation flow domains --------------------

2
Cantera/clib/src/ctonedim.h
View File

@@ -27,11 +27,13 @@ extern "C" {
double DLL_IMPORT bdry_temperature(int i);
double DLL_IMPORT bdry_massFraction(int i, int k);
double DLL_IMPORT bdry_mdot(int i);
int DLL_IMPORT reactingsurf_setkineticsmgr(int i, int j);
int DLL_IMPORT inlet_new();
int DLL_IMPORT outlet_new();
int DLL_IMPORT symm_new();
int DLL_IMPORT surf_new();
int DLL_IMPORT reactingsurf_new();
int DLL_IMPORT stflow_new(int iph, int ikin, int itr);
int DLL_IMPORT stflow_setPressure(int i, double p);

4
Cantera/matlab/cantera/examples/flame.m
View File

@@ -64,8 +64,8 @@ end
zz = z(flow);
dz = zz(end) - zz(1);
setProfile(f, 2, {'u', 'V'}, [0.0 1.0
mdot0/rho0 -mdot1/rho0
0.0 0.0]);
mdot0/rho0 -mdot1/rho0
0.0 0.0]);
setProfile(f, 2, 'T', [0.0 z1 1.0; t0 2000.0 t1]);
for n = 1:nSpecies(gas)

7
Cantera/matlab/cantera/private/onedimmethods.cpp
View File

@@ -57,6 +57,13 @@ void onedimmethods( int nlhs, mxArray *plhs[],
indx = outlet_new();
break;
// construct a new ReactingSurf1D instance
case 6:
checkNArgs(4, nrhs);
indx = reactingsurf_new();
reactingsurf_setkineticsmgr(indx, getInt(prhs[3]));
break;
// construct a new Sim1D instance
case 8:
checkNArgs(5, nrhs);

87
Cantera/python/Cantera/ctml_writer.py
View File

@@ -49,6 +49,15 @@ _reactions = []
_atw = {}
_mw = {}
_valsp = ''
_valrxn = ''
def validate(species = 'yes', reactions = 'yes'):
global _valsp
global _valrxn
_valsp = species
_valrxn = reactions
def isnum(a):
"""True if a is an integer or floating-point number."""
if type(a) == types.IntType or type(a) == types.FloatType:
@@ -119,6 +128,10 @@ def ufmt(base, n):
def write():
"""write the CTML file."""
x = XML_Node("ctml")
v = x.addChild("validate")
v["species"] = _valsp
v["reactions"] = _valrxn
for ph in _phases:
ph.build(x)
s = species_set(name = _name, species = _species)
@@ -264,6 +277,8 @@ class species(writer):
global _species
_species.append(self)
global _speciesnames
if name in _speciesnames:
raise CanteraError('species '+name+' multiply defined.')
_speciesnames.append(name)
@@ -437,11 +452,16 @@ class reaction(writer):
def __init__(self,
equation = '',
kf = None,
id = '',
id = '',
special = []
):
self._id = id
self._e = equation
if type(special) == types.StringType:
self._special = [special]
else:
self._special = special
global _reactions
self._num = len(_reactions)+1
r = ''
@@ -493,37 +513,10 @@ class reaction(writer):
if ph.is_ideal_gas():
self._igspecies.append(s)
break
## if nm < 0:
## if _handle_undeclared_species == _WARN:
## print 'Warning... skipping reaction '+self._e
## print 'because species '+s+' is undeclared'
## if _handle_undeclared_species <= _WARN:
## return 0
## else:
## raise CanteraError('undeclared species '+s+' in reaction '+self._e)
## else:
mdim += nm*ns
ldim += nl*ns
## for s in self._p.keys():
## ns = self._p[s]
## ok = 0
## for ph in _phases:
## if ph.has_species(s):
## ok = 1
## break
## if not ok:
## if _handle_undeclared_species == _WARN:
## print 'Warning... skipping reaction '+self._e
## print 'because species '+s+' is undeclared'
## if _handle_undeclared_species <= _WARN:
## return 0
## else:
## raise CanteraError('undeclared species '+s+' in reaction '+self._e)
p.addComment(" reaction "+id+" ")
r = p.addChild('reaction')
r['id'] = id
@@ -531,6 +524,13 @@ class reaction(writer):
r['reversible'] = 'yes'
else:
r['reversible'] = 'no'
nspecial = len(self._special)
for nss in range(nspecial):
s = self._special[nss]
if s == 'duplicate':
r['duplicate'] = 'yes'
elif s == 'negative_A':
r['negative_A'] = 'yes'
ee = self._e.replace('<','[')
ee = ee.replace('>',']')
@@ -596,15 +596,17 @@ class reaction(writer):
#-------------------
class three_body_reaction(reaction):
def __init__(self,
equation = '',
kf = None,
efficiencies = '',
id = '',
id = '',
special = []
):
reaction.__init__(self, equation, kf, id)
reaction.__init__(self, equation, kf, id, special)
self._type = 'threeBody'
self._effm = 1.0
self._eff = efficiencies
@@ -615,8 +617,7 @@ class three_body_reaction(reaction):
del self._r[r]
for p in self._p.keys():
if p == 'M' or p == 'm':
del self._p[p]
del self._p[p]
def build(self, p):
@@ -639,7 +640,8 @@ class falloff_reaction(reaction):
kf = None,
efficiencies = '',
falloff = None,
id = ''
id = '',
special = []
):
kf2 = (kf, kf0)
@@ -693,7 +695,8 @@ class surface_reaction(reaction):
equation = '',
kf = None,
stick = None,
id = ''):
id = '',
special = []):
if stick:
reaction.__init__(self, equation, stick, id)
else:
@@ -775,26 +778,29 @@ class phase(writer):
sptoks = spnames.replace(',',' ').split()
for s in sptoks:
if self._spmap.has_key(s):
raise 'multiply-defined species '+s+' in phase '+self._name
raise CanteraError('Multiply-declared species '+s+' in phase '+self._name)
self._spmap[s] = self._dim
self._rxns = reactions
# check that species have been declared
if len(self._spmap) == 0:
raise 'No species declared for phase '+self._name
raise CanteraError('No species declared for phase '+self._name)
# and that only one species is declared if it is a pure phase
if self.is_pure() and len(self._spmap) > 1:
raise 'Pure phases may only declare one species, but phase '+self._name+' declares '+`len(self._spmap)`+'.'
raise CanteraError('Stoichiometric phases must declare exactly one species, \n'+
'but phase '+self._name+' declares '+`len(self._spmap)`+'.')
self._initial = initial_state
# add this phase to the global phase list
global _phases
_phases.append(self)
def is_ideal_gas(self):
"""True if the entry represents an ideal gas."""
return 0
def is_pure(self):
@@ -853,6 +859,12 @@ class phase(writer):
ph = p.addChild('phase')
ph['id'] = self._name
ph['dim'] = `self._dim`
# ------- error tests -------
err = ph.addChild('validation')
err.addChild('duplicateReactions','halt')
err.addChild('thermo','warn')
e = ph.addChild('elementArray',self._el)
e['datasrc'] = 'elements.xml'
for s in self._sp:
@@ -1116,3 +1128,4 @@ class Lindemann:
get_atomic_wts()
validate()

3
Cantera/python/Cantera/solution.py
View File

@@ -39,8 +39,7 @@ class Solution(ThermoPhase, Kinetics, Transport):
id = ""
if len(fn) > 1:
id = fn[1]
fn = fn[0]
fn = fn[0]
fname = os.path.basename(fn)
ff = os.path.splitext(fname)

167
Cantera/src/importCTML.cpp
View File

@@ -50,6 +50,12 @@ using namespace ctml;
GasKineticsWriter* writer = 0;
vector< map<int, doublereal> > _reactiondata;
vector<string> _eqn;
vector_int _dup;
vector<bool> _rev;
namespace Cantera {
/*
* First we define a coule of typedef's which will
@@ -69,7 +75,7 @@ namespace Cantera {
* Install a NASA polynomial thermodynamic property
* parameterization for species k.
*/
void installNasaThermo(SpeciesThermo& sp, int k, XML_Node& f0,
static void installNasaThermo(SpeciesThermo& sp, int k, XML_Node& f0,
XML_Node& f1) {
doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
@@ -113,7 +119,7 @@ namespace Cantera {
* Install a Shomate polynomial thermodynamic property
* parameterization for species k.
*/
void installShomateThermo(SpeciesThermo& sp, int k, XML_Node& f) {
static void installShomateThermo(SpeciesThermo& sp, int k, XML_Node& f) {
doublereal tmin, tmid, tmax;
tmin = fpValue(f["Tmin"]);
tmid = fpValue(f["Tmid"]);
@@ -139,17 +145,15 @@ namespace Cantera {
}
/**
* Install a Shomate polynomial thermodynamic property
* Install a constant-cp thermodynamic property
* parameterization for species k.
*/
void installSimpleThermo(SpeciesThermo& sp, int k, XML_Node& f) {
static void installSimpleThermo(SpeciesThermo& sp, int k, XML_Node& f) {
doublereal tmin, tmax;
tmin = fpValue(f["Tmin"]);
tmax = fpValue(f["Tmax"]);
if (tmax == 0.0) tmax = 1.0e30;
//map<string, doublereal> fp;
//getFloats(f, fp);
vector_fp c(4);
c[0] = getFloat(f, "t0", "-");
c[1] = getFloat(f, "h0", "-");
@@ -163,7 +167,7 @@ namespace Cantera {
* Install a species into a ThermoPhase object, which defines
* the phase thermodynamics and speciation
*/
bool installSpecies(int k, XML_Node& s, thermo_t& p,
static bool installSpecies(int k, XML_Node& s, thermo_t& p,
SpeciesThermo& spthermo, int rule) {
// get the composition of the species
@@ -179,9 +183,8 @@ namespace Cantera {
if (rule == 0)
throw
CanteraError("installSpecies",
"species " + s["name"] +
" contains undeclared element " +
_b->first);
"Species " + s["name"] +
" contains undeclared element " + _b->first);
else
return false;
}
@@ -202,8 +205,7 @@ namespace Cantera {
doublereal sz = 1.0;
if (s.hasChild("size")) sz = getFloat(s, "size");
p.addUniqueSpecies(s["name"], ecomp.begin(),
chrg, sz);
p.addUniqueSpecies(s["name"], ecomp.begin(), chrg, sz);
// get thermo
XML_Node& thermo = s.child("thermo");
@@ -211,9 +213,6 @@ namespace Cantera {
int nc = tp.size();
if (nc == 1) {
XML_Node& f = *tp[0];
//if (f.name() == "NASA") {
// installNasaThermo(spthermo, k, f);
//}
if (f.name() == "Shomate") {
installShomateThermo(spthermo, k, f);
}
@@ -268,7 +267,7 @@ namespace Cantera {
* rule = If we fail to find a species, we will throw an error
* if rule != 1.
*/
bool getReagents(XML_Node& rxn, kinetics_t& kin, int rp,
static bool getReagents(XML_Node& rxn, kinetics_t& kin, int rp,
string default_phase,
vector_int& spnum, vector_int& stoich, vector_fp& order,
int rule) {
@@ -312,7 +311,7 @@ namespace Cantera {
return false;
else {
throw CanteraError("getReagents",
"undeclared reactant or product species "+sp);
"Undeclared reactant or product species "+sp);
return false;
}
}
@@ -337,7 +336,7 @@ namespace Cantera {
* Arrhenius expression is
* k = A T^(b) exp (-Ea / RT).
*/
void getArrhenius(XML_Node& node, int& highlow, doublereal& A,
static void getArrhenius(XML_Node& node, int& highlow, doublereal& A,
doublereal& b, doublereal& E) {
if (node["name"] == "k0")
@@ -366,9 +365,7 @@ namespace Cantera {
const ThermoPhase& p = kin.thermo(np);
klocal = p.speciesIndex(kin.kineticsSpeciesName(k));
if (p.eosType() == cSurf) {
cout << "f before = " << f << endl;
f /= pow(p.standardConcentration(klocal),ns);
cout << "f, k, klocal = " << f << " " << k << " " << klocal << endl;
}
else
not_surf++;
@@ -383,14 +380,13 @@ namespace Cantera {
b = getFloat(node, "b") + 0.5;
E = getFloat(node, "E", "actEnergy");
E /= GasConstant;
cout << "A, cbar, f, b, E: " << A << " " << cbar << " " << f << " " << b << " " << E << endl;
}
/**
* Get falloff parameters for a reaction.
*/
void getFalloff(node_t& f, ReactionData& rdata) {
static void getFalloff(node_t& f, ReactionData& rdata) {
string type = f["type"];
vector<string> p;
getStringArray(f,p);
@@ -415,7 +411,7 @@ namespace Cantera {
* reaction mechanism is homogeneous, so that all species belong
* to phase(0) of 'kin'.
*/
void getEfficiencies(node_t& eff, kinetics_t& kin, ReactionData& rdata) {
static void getEfficiencies(node_t& eff, kinetics_t& kin, ReactionData& rdata) {
// set the default collision efficiency
rdata.default_3b_eff = fpValue(eff["default"]);
@@ -440,8 +436,8 @@ namespace Cantera {
* This function will fill in more fields in the ReactionData object.
*
*/
void getRateCoefficient(node_t& kf, kinetics_t& kin,
ReactionData& rdata) {
static void getRateCoefficient(node_t& kf, kinetics_t& kin,
ReactionData& rdata, int negA) {
int nc = kf.nChildren();
const nodeset_t& kf_children = kf.children();
@@ -459,6 +455,10 @@ namespace Cantera {
|| rdata.reactionType == ELEMENTARY_RXN)
chigh = coeff;
else clow = coeff;
if (coeff[0] <= 0.0 && negA == 0) {
throw CanteraError("getRateCoefficient",
"negative or zero A coefficient for reaction "+int2str(rdata.number));
}
}
else if (nm == "Stick") {
vector_fp coeff(3);
@@ -467,6 +467,10 @@ namespace Cantera {
int isp = th.speciesIndex(spname);
double mw = th.molecularWeights()[isp];
getStick(c, mw, kin, rdata, coeff[0], coeff[1], coeff[2]);
if (coeff[0] <= 0.0 && negA == 0) {
throw CanteraError("getRateCoefficient",
"negative or zero A coefficient for reaction "+int2str(rdata.number));
}
chigh = coeff;
}
@@ -761,6 +765,29 @@ namespace Cantera {
return true;
}
/**
* This function returns true if two reactions are duplicates of
* one another, and false otherwise. The input arguments are two
* maps from species number to stoichiometric coefficient, one for
* each reaction. The reactions are considered duplicates if their
* stoichiometric coefficients have the same ratio for all
* species.
*/
static bool isDuplicateReaction(map<int, doublereal>& r1, map<int, doublereal>& r2) {
map<int, doublereal>::const_iterator b = r1.begin(), e = r1.end();
int k1 = b->first;
doublereal ratio = r2[k1]/r1[k1];
if (r1[k1] == 0.0 || r2[k1] == 0.0) return false;
++b;
for (; b != e; ++b) {
k1 = b->first;
if (r1[k1] == 0.0 || r2[k1] == 0.0) return false;
if (fabs(r2[k1]/r1[k1] - ratio) > 1.e-8) return false;
}
return true;
}
/**
* Install an individual reaction into the kinetics mechanism
@@ -773,8 +800,8 @@ namespace Cantera {
* rule = Provides a rule for specifying how to handle reactions
* which involve missing species.
*/
bool installReaction(int i, XML_Node& r, Kinetics* k,
string default_phase, int rule) {
static bool installReaction(int i, XML_Node& r, Kinetics* k,
string default_phase, int rule, bool check_for_duplicates) {
Kinetics& kin = *k;
/*
@@ -790,11 +817,18 @@ namespace Cantera {
int nn, eqlen;
vector_fp dummy;
int dup = 0;
if (r.hasAttrib("duplicate")) dup = 1;
int negA = 0;
if (r.hasAttrib("negative_A")) negA = 1;
/*
* This seemingly simple expression goes and finds the child element,
* "equation". Then it treats all of the contents of the "equation"
* as a string, and returns it the variable eqn. We post process
* the string to get rid of [ and ] characters for some reason.
* the string to convert [ and ] characters into < and >, which cannot be
* stored in an XML file.
*/
if (r.hasChild("equation"))
eqn = r("equation");
@@ -821,6 +855,47 @@ namespace Cantera {
return false;
}
string isrev = r["reversible"];
if (isrev == "yes" || isrev == "true")
rdata.reversible = true;
if (check_for_duplicates) {
doublereal c = 0.0;
map<int, doublereal> rxnstoich;
rxnstoich.clear();
int nr = rdata.reactants.size();
for (nn = 0; nn < nr; nn++) {
rxnstoich[rdata.reactants[nn]] -= rdata.rstoich[nn];
}
int np = rdata.products.size();
for (nn = 0; nn < np; nn++) {
rxnstoich[rdata.products[nn]] += rdata.pstoich[nn];
}
int nrxns = _reactiondata.size();
for (nn = 0; nn < nrxns; nn++) {
c = isDuplicateReaction(rxnstoich, _reactiondata[nn]);
if (c > 0.0
|| (c < 0.0 && rdata.reversible)
|| (c < 0.0 && _rev[nn])) {
if ((!dup || !_dup[nn])) {
string msg = string("Undeclared duplicate reactions detected: \n")
+"Reaction "+int2str(nn+1)+": "+_eqn[nn]
+"\nReaction "+int2str(i+1)+": "+eqn+"\n";
_reactiondata.clear();
_eqn.clear();
_rev.clear();
_dup.clear();
throw CanteraError("installReaction",msg);
}
}
}
_dup.push_back(dup);
_rev.push_back(rdata.reversible);
_eqn.push_back(eqn);
_reactiondata.push_back(rxnstoich);
}
rdata.equation = eqn;
rdata.reversible = false;
rdata.number = i;
@@ -849,11 +924,7 @@ namespace Cantera {
throw CanteraError("installReaction",
"Unknown reaction type: " + typ);
string isrev = r["reversible"];
if (isrev == "yes" || isrev == "true")
rdata.reversible = true;
getRateCoefficient(r.child("rateCoeff"), kin, rdata);
getRateCoefficient(r.child("rateCoeff"), kin, rdata, negA);
/*
* Ok we have read everything in about the reaction. Add it
* to the kinetics object by calling the Kinetics member function,
@@ -878,7 +949,13 @@ namespace Cantera {
* If there is a problem, return false.
*/
bool installReactionArrays(XML_Node& p, Kinetics& kin,
string default_phase) {
string default_phase, bool check_for_duplicates) {
_eqn.clear();
_dup.clear();
_reactiondata.clear();
_rev.clear();
vector<XML_Node*> rarrays;
int itot = 0;
/*
@@ -944,17 +1021,13 @@ namespace Cantera {
XML_Node* r = allrxns[i];
if (r) {
if (installReaction(itot, *r, &kin,
default_phase, rxnrule)) ++itot;
default_phase, rxnrule, check_for_duplicates)) ++itot;
}
}
}
else {
for (int nii = 0; nii < ninc; nii++) {
XML_Node& ii = *incl[nii];
//vector<string> rxn_ids;
//string pref = ii["prefix"];
//int imin = atoi(ii["min"].c_str());
//int imax = atoi(ii["max"].c_str());
string imin = ii["min"];
string imax = ii["max"];
for (i = 0; i < nrxns; i++) {
@@ -962,8 +1035,6 @@ namespace Cantera {
string rxid;
if (r) {
rxid = (*r)["id"];
//cout << rxid << " " << imin << " " << imax << endl;
//cout << (rxid >= imin) << " " << (rxid <= imax) << endl;
/*
* To decide whether the reaction is included or not
* we do a lexical min max and operation. This
@@ -971,19 +1042,23 @@ namespace Cantera {
*/
if ((rxid >= imin) && (rxid <= imax)) {
if (installReaction(itot, *r, &kin,
default_phase, rxnrule)) ++itot;
default_phase, rxnrule, check_for_duplicates)) ++itot;
}
}
}
}
}
}
/*
* Finalize the installation of the kinetics, now that we know
* the true number of reactions in the mechanism, itot.
*/
kin.finalize();
writer = 0;
_eqn.clear();
_dup.clear();
_reactiondata.clear();
return true;
}
@@ -999,6 +1074,12 @@ namespace Cantera {
// This phase will be the default one
string default_phase = phase["id"];
bool check_for_duplicates = false;
if (phase.parent()->hasChild("validate")) {
XML_Node& d = phase.parent()->child("validate");
if (d["reactions"] == "yes") check_for_duplicates = true;
}
// if other phases are involved in the reaction mechanism,
// they must be listed in a 'phaseArray' child
// element. Homogeneous mechanisms do not need to include a
@@ -1049,7 +1130,7 @@ namespace Cantera {
kin.init();
// Install the reactions.
return installReactionArrays(phase, kin, default_phase);
return installReactionArrays(phase, kin, default_phase, check_for_duplicates);
}
/**

2
Cantera/src/importCTML.h
View File

@@ -25,7 +25,7 @@ namespace Cantera {
bool importKinetics(XML_Node& phase, vector<thermophase_t*> th,
Kinetics* kin);
bool installReactionArrays(XML_Node& parent, Kinetics& kin,
string default_phase);
string default_phase, bool check_for_duplicates = false);
ThermoPhase* newPhase(XML_Node& phase);
bool buildSolutionFromXML(XML_Node& root, string id, string nm,
ThermoPhase* th, Kinetics* k);

73
Cantera/src/oneD/Inlet1D.h
View File

@@ -109,7 +109,7 @@ namespace Cantera {
* without parameters, a left inlet (facing right) is
* constructed).
*/
Inlet1D() : m_V0(0.0), m_nsp(0), m_flow(0) {
Inlet1D() : Bdry1D(), m_V0(0.0), m_nsp(0), m_flow(0) {
m_type = cInletType;
m_xstr = "";
}
@@ -200,7 +200,7 @@ namespace Cantera {
public:
Symm1D() {
Symm1D() : Bdry1D() {
m_type = cSymmType;
}
virtual ~Symm1D(){}
@@ -228,7 +228,7 @@ namespace Cantera {
public:
Outlet1D() {
Outlet1D() : Bdry1D() {
m_type = cOutletType;
}
virtual ~Outlet1D(){}
@@ -260,7 +260,7 @@ namespace Cantera {
public:
Surf1D() {
Surf1D() : Bdry1D() {
m_type = cSurfType;
}
virtual ~Surf1D(){}
@@ -299,6 +299,71 @@ namespace Cantera {
};
/**
* A reacting surface.
*
*/
class ReactingSurf1D : public Bdry1D {
public:
ReactingSurf1D() : Bdry1D(),
m_kin(0), m_surfindex(0), m_nsp(0) {
m_type = cSurfType;
}
void setKineticsMgr(InterfaceKinetics* kin) {
m_surfindex = kin->surfacePhaseIndex();
m_sphase = (SurfPhase*)&kin->thermo(m_surfindex);
m_nsp = m_sphase->nSpecies();
m_enabled = true;
}
virtual ~ReactingSurf1D(){}
virtual string componentName(int n) const;
virtual void init();
virtual void eval(int jg, doublereal* xg, doublereal* rg,
integer* diagg, doublereal rdt);
virtual void save(XML_Node& o, doublereal* soln);
virtual void restore(XML_Node& dom, doublereal* soln);
virtual void _getInitialSoln(doublereal* x) {
x[0] = m_temp;
m_sphase->getCoverages(x+1);
}
virtual void _finalize(const doublereal* x) {
; //m_temp = x[0];
}
virtual void showSolution(const doublereal* x) {
char buf[80];
sprintf(buf, " Temperature: %10.4g K \n", x[0]);
writelog(buf);
writelog(" Coverages: \n");
for (int k = 0; k < m_nsp; k++) {
sprintf(buf, " %20s %10.4g \n", m_sphase->speciesName(k).c_str(),
x[k+1]);
writelog(buf);
}
writelog("\n");
}
protected:
InterfaceKinetics* m_kin;
SurfPhase* m_sphase;
int m_surfindex, m_nsp;
bool m_enabled;
vector_fp m_fixed_cov;
vector_fp m_work;
};
}
// // };

2
Cantera/src/oneD/Sim1D.cpp
View File

@@ -337,7 +337,7 @@ namespace Cantera {
}
}
else {
throw CanteraError("refine","keepPoint is false at m = "+int2str(m));
; // throw CanteraError("refine","keepPoint is false at m = "+int2str(m));
}
}
dsize.push_back(znew.size() - nstart);

116
Cantera/src/oneD/Solid1D.cpp
View File

@@ -31,32 +31,23 @@ namespace Cantera {
int Solid1D::c_T_loc = 0;
int Solid1D::c_phi_loc = 1;
int Solid1D::c_Y_loc = 2;
int Solid1D::c_C_loc = 1;
Solid1D::Solid1D(ThermoPhase* ph, int nsp, int points) :
Domain1D(nsp+1, points),
m_nsp(nsp),
m_thermo(0),
Solid1D::Solid1D(ThermoPhase* ph, int points) :
Domain1D(1, points),
m_kin(0),
m_trans(0),
m_jac(0),
m_ok(false),
m_ok(false)
{
m_type = cSolidType;
m_points = points;
m_thermo = ph;
if (ph == 0) return; // used to create a dummy object
int nsp2 = m_thermo->nSpecies();
if (nsp2 != m_nsp) {
m_nsp = nsp2;
Domain1D::resize(m_nsp+1, points);
}
if (ph == 0) { m_nsp = 1; return; }// used to create a dummy object
m_nsp = m_thermo->nSpecies();
Domain1D::resize(m_nsp+1, points);
// make a local copy of the species molecular weight vector
m_wt = m_thermo->molecularWeights();
@@ -65,13 +56,12 @@ namespace Cantera {
// turn off the energy equation at all points
m_do_energy.resize(m_points,false);
m_do_gauss.resize(m_points,false);
m_do_species.resize(m_nsp,false);
m_diff.resize(m_nsp*m_points);
m_flux.resize(m_nsp,m_points);
m_wdot.resize(m_nsp,m_points, 0.0);
m_ybar.resize(m_nsp);
m_cbar.resize(m_nsp);
//-------------- default solution bounds --------------------
@@ -82,11 +72,11 @@ namespace Cantera {
vmin[c_T_loc] = 200.0;
vmax[c_T_loc]= 1.e9;
// mass fraction bounds
// concentration bounds
int k;
for (k = 0; k < m_nsp; k++) {
vmin[c_Y_loc + k] = -1.0e-5;
vmax[c_Y_loc + k] = 1.0e5;
vmin[c_C_loc + k] = -1.0e-5;
vmax[c_C_loc + k] = 1.0e5;
}
setBounds(vmin.size(), vmin.begin(), vmax.size(), vmax.begin());
@@ -117,18 +107,13 @@ namespace Cantera {
m_wtm.resize(m_points, 0.0);
m_cp.resize(m_points, 0.0);
m_tcon.resize(m_points, 0.0);
m_cdens.resize(m_points, 0.0);
m_diff.resize(m_nsp*m_points);
m_flux.resize(m_nsp,m_points);
m_wdot.resize(m_nsp,m_points, 0.0);
m_do_energy.resize(m_points,false);
m_do_gauss.resize(m_points,false);
m_do_species.resize(m_nsp,false);
m_fixedtemp.resize(m_points);
m_fixedphi.resize(m_points);
m_dz.resize(m_points-1);
m_z.resize(m_points);
@@ -164,9 +149,8 @@ namespace Cantera {
*/
void Solid1D::setThermoState(const doublereal* x,int j) {
m_thermo->setTemperature(T(x,j));
m_thermo->setElectricPotential(phi(x,j));
const doublereal* yy = x + m_nv*j + 1;
m_thermo->setMassFractions_NoNorm(yy);
m_thermo->setConcentrations(yy);
}
@@ -176,12 +160,11 @@ namespace Cantera {
*/
void Solid1D::setStateAtMidpoint(const doublereal* x,int j) {
m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
m_thermo->setElectricPotential(0.5*(phi(x,j)+phi(x,j+1)));
const doublereal* yyj = x + m_nv*j + 1;
const doublereal* yyjp = x + m_nv*(j+1) + 1;
const doublereal* ccj = x + m_nv*j + 1;
const doublereal* ccjp = x + m_nv*(j+1) + 1;
for (int k = 0; k < m_nsp; k++)
m_ybar[k] = 0.5*(yyj[k] + yyjp[k]);
m_thermo->setMassFractions_NoNorm(m_ybar.begin());
m_ybar[k] = 0.5*(ccj[k] + ccjp[k]);
m_thermo->setConcentrations(m_cbar.begin());
}
@@ -221,7 +204,6 @@ namespace Cantera {
int j, k;
//-----------------------------------------------------
// update properties
//-----------------------------------------------------
@@ -242,7 +224,6 @@ namespace Cantera {
for (j = j0; j <= j1; j++) {
setThermoState(j);
m_cdens[j] = m_thermo->chargeDensity();
}
//----------------------------------------------------
@@ -259,13 +240,12 @@ namespace Cantera {
if (j == 0) {
rsd[index(c_T_loc,0)] = T(x,0);
rsd[index(c_phi_loc,0)] = phi(x,0);
// The default boundary condition for species is zero
// flux. However, the boundary object may modify
// this.
for (k = 0; k < m_nsp; k++) {
rsd[index(c_Y_loc + k, 0)] = - m_flux(k,0);
rsd[index(c_C_loc + k, 0)] = - m_flux(k,0);
}
}
@@ -278,14 +258,9 @@ namespace Cantera {
else if (j == m_points - 1) {
rsd[index(c_T_loc,j)] = T(x,j);
rsd[index(c_phi_loc,j)] = phi(x,j);
doublereal sum = 0.0;
for (k = 0; k < m_nsp; k++) {
sum += Y(x,k,j);
rsd[index(k+c_Y_loc,j)] = m_flux(k,j-1);
rsd[index(k+c_C_loc,j)] = m_flux(k,j-1);
}
rsd[index(c_Y_loc,j)] = 1.0 - sum;
diag[index(c_Y_loc,j)] = 0;
}
@@ -303,15 +278,14 @@ namespace Cantera {
//-------------------------------------------------
getWdot(x,j);
doublereal convec, diffus;
doublereal diffus;
for (k = 0; k < m_nsp; k++) {
diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))
/(z(j+1) - z(j-1));
rsd[index(c_Y_loc + k, j)]
= (m_wt[k]*(wdot(k,j) )
- diffus)/m_rho[j]
- rdt*(Y(x,k,j) - Y_prev(k,j));
diag[index(c_Y_loc + k, j)] = 1;
rsd[index(c_C_loc + k, j)]
= wdot(k,j) - diffus
- rdt*(C(x,k,j) - C_prev(k,j));
diag[index(c_C_loc + k, j)] = 1;
}
//-----------------------------------------------
@@ -326,12 +300,6 @@ namespace Cantera {
rsd[index(c_T_loc, j)] -= rdt*(T(x,j) - T_prev(j));
diag[index(c_T_loc, j)] = 1;
}
//----------------------------------------------
// Gauss's equation
//----------------------------------------------
rsd[index(c_phi_loc, j)] = m_cdens[j] - divDisplCurr(x,j);
}
// residual equations if the energy or species equations
@@ -341,10 +309,6 @@ namespace Cantera {
rsd[index(c_T_loc, j)] = T(x,j) - T_fixed(j);
diag[index(c_T_loc, j)] = 0;
}
if (!m_do_gauss[j]) {
rsd[index(c_phi_loc, j)] = phi(x,j) - phi_fixed(j);
diag[index(c_phi_loc, j)] = 0;
}
}
}
@@ -424,19 +388,14 @@ namespace Cantera {
doublereal dphidz, a1;
for (j = j0; j < j1; j++) {
sum = 0.0;
wtm = m_wtm[j];
rho = density(j);
dz = z(j+1) - z(j);
dphidz = (phi(x,j+1) - phi(x,j))/dz;
a1 = rho*Faraday*dphidz/(GasConstant * T(x,j));
for (k = 0; k < m_nsp; k++) {
m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
m_flux(k,j) += a1*0.5*(Y(x,k,j)
+ Y(x,k,j+1))*m_diff[k+m_nsp*j]*m_charge[k];
m_flux(k,j) = m_diff[k+m_nsp*j] *
(C(x,k,j) - C(x,k,j+1))/dz;
sum -= m_flux(k,j);
}
for (k = 0; k < m_nsp; k++) m_flux(k,j) += Y(x,k,j)*sum;
for (k = 0; k < m_nsp; k++) m_flux(k,j) += C(x,k,j)*sum;
}
break;
}
@@ -448,14 +407,13 @@ namespace Cantera {
s << "TITLE = \"" + title + "\"" << endl;
s << "VARIABLES = \"Z (m)\"" << endl;
s << "\"T (K)\"" << endl;
s << "\"phi (V)\"" << endl;
for (k = 0; k < m_nsp; k++) {
s << "\"" << m_thermo->speciesName(k) << "\"" << endl;
}
s << "ZONE T=\"c" << zone << "\"" << endl;
s << " I=" << m_points << ",J=1,K=1,F=POINT" << endl;
s << "DT=(SINGLE SINGLE";
s << "DT=(SINGLE";
for (k = 0; k < m_nsp; k++) s << " SINGLE";
s << " )" << endl;
for (j = 0; j < m_points; j++) {
@@ -471,9 +429,8 @@ namespace Cantera {
string Solid1D::componentName(int n) const {
switch(n) {
case c_T_loc: return "T";
case c_phi_loc: return "phi";
default:
if (n >= (int) 1 && n < (int) (c_Y_loc + m_nsp)) {
if (n >= (int) 1 && n < (int) (c_C_loc + m_nsp)) {
return m_thermo->speciesName(n - 1);
}
else
@@ -548,21 +505,13 @@ namespace Cantera {
}
else goto error;
}
else if (nm == "phi") {
writelog("potential ");
if ((int) x.size() == np) {
for (j = 0; j < np; j++)
soln[index(c_phi_loc,j)] = x[j];
}
else goto error;
}
else if (m_thermo->speciesIndex(nm) >= 0) {
writelog(nm+" ");
if ((int) x.size() == np) {
k = m_thermo->speciesIndex(nm);
did_species[k] = 1;
for (j = 0; j < np; j++)
soln[index(k+c_Y_loc,j)] = x[j];
soln[index(k+c_C_loc,j)] = x[j];
}
}
else
@@ -615,13 +564,10 @@ namespace Cantera {
soln.getRow(c_T_loc,x.begin());
addFloatArray(gv,"T",x.size(),x.begin(),"K","temperature",0.0);
soln.getRow(c_phi_loc,x.begin());
addFloatArray(gv,"phi",x.size(),x.begin(),"V","potential",0.0);
for (k = 0; k < m_nsp; k++) {
soln.getRow(c_Y_loc+k,x.begin());
soln.getRow(c_C_loc+k,x.begin());
addFloatArray(gv,m_thermo->speciesName(k),
x.size(),x.begin(),"","massFraction",0.0,1.0);
x.size(),x.begin(),"","concentration",0.0,1.0);
}
}

195
Cantera/src/oneD/boundaries1D.cpp
View File

@@ -469,10 +469,203 @@ namespace Cantera {
soln[0] = x["temperature"];
resize(1,1);
}
}
//-----------------------------------------------------------
//
// ReactingSurf1D
//
//-----------------------------------------------------------
string ReactingSurf1D::componentName(int n) const {
if (n == 0) return "temperature";
else if (n < m_nsp + 1)
return m_sphase->speciesName(n-1);
else
return "<unknown>";
}
void ReactingSurf1D::
init() {
_init(m_nsp+1);
m_fixed_cov.resize(m_nsp, 0.0);
m_fixed_cov[0] = 1.0;
m_work.resize(m_kin->nTotalSpecies());
// 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;
}
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;
atol[n] = 1.0e-9;
}
atol[0] = 1.0e-4;
setTolerances(m_nv, rtol.begin(), m_nv, atol.begin());
}
void ReactingSurf1D::
eval(int jg, doublereal* xg, doublereal* rg,
integer* diagg, doublereal rdt) {
if (jg >= 0 && (jg < firstPoint() - 2 || jg > lastPoint() + 2)) return;
// start of local part of global arrays
doublereal* x = xg + loc();
doublereal* r = rg + loc();
integer* diag = diagg + loc();
doublereal *xb, *rb;
// specified surface temp
r[0] = x[0] - m_temp;
// set the coverages
doublereal sum = 0.0;
int k;
for (k = 0; k < m_nsp; k++) {
m_work[k] = x[k];
sum += x[k];
}
m_sphase->setCoverages(m_work.begin());
// set the left gas state to the adjacent point
int leftloc = 0, rightloc = 0;
int pnt = 0;
if (m_flow_left) {
leftloc = m_flow_left->loc();
pnt = m_flow_left->nPoints() - 1;
m_flow_left->setGas(xg + leftloc, pnt);
}
if (m_flow_right) {
rightloc = m_flow_right->loc();
m_flow_right->setGas(xg + rightloc, 0);
}
m_kin->getNetProductionRates(m_work.begin());
doublereal rs0 = 1.0/m_sphase->siteDensity();
//scale(m_work.begin(), m_work.end(), m_work.begin(), m_mult[0]);
bool enabled = true;
int ioffset = m_kin->kineticsSpeciesIndex(0, m_surfindex);
if (m_enabled) {
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[imx] = 1.0 - sum;
diag[imx] = 0;
}
else {
r[k] = x[k] - m_fixed_cov[k];
diag[k] = 0;
}
// gas-phase residuals
// doublereal rho;
// if (m_flow_left) {
// rho = m_phase_left->density();
// doublereal rdz = 2.0/
// (m_flow_left->z(m_left_points-1) -
// m_flow_left->z(m_left_points - 2));
// for (k = 0; k < m_left_nsp; k++)
// m_work[k + m_start_left] *= m_molwt_left[k];
// int ileft = loc() - m_left_nv;
// // if the energy equation is enabled at this point,
// // set the gas temperature to the surface temperature
// if (m_flow_left->doEnergy(pnt)) {
// rg[ileft + 2] = xg[ileft + 2] - m_sphase->temperature();
// }
// for (k = 1; k < m_left_nsp; k++) {
// if (enabled && m_flow_left->doSpecies(k)) {
// rg[ileft + 4 + k] += m_work[k + m_start_left];
// //+= rdz*m_work[k + m_sp_left]/rho;
// }
// }
// }
// if (m_flow_right) {
// for (k = 0; k < m_right_nsp; k++)
// m_work[k + m_start_right] *= m_molwt_right[k];
// int iright = loc() + m_nsp;
// rg[iright + 2] -= m_sphase->temperature();
// //r[iright + 3] = x[iright];
// for (k = 0; k < m_right_nsp; k++) {
// rg[iright + 4 + k] -= m_work[k + m_start_right];
// }
// }
// diag[0] = 0;
// int nc;
// if (m_flow_right) {
// rb = r + 1;
// xb = x + 1;
// rb[2] = xb[2] - x[0]; // specified T
// }
// if (m_flow_left) {
// nc = m_flow_left->nComponents();
// rb = r - nc;
// xb = x - nc;
// rb[2] = xb[2] - x[0]; // specified T
// }
}
void ReactingSurf1D::
save(XML_Node& o, doublereal* soln) {
doublereal* s = soln + loc();
//XML_Node& inlt = o.addChild("inlet");
XML_Node& inlt = o.addChild("domain");
inlt.addAttribute("id",id());
inlt.addAttribute("points",1);
inlt.addAttribute("type","surface");
inlt.addAttribute("components",nComponents());
for (int k = 0; k < nComponents(); k++) {
ctml::addFloat(inlt, componentName(k), s[k], "", "",0.0, 1.0);
}
}
void ReactingSurf1D::
restore(XML_Node& dom, doublereal* soln) {
map<string, double> x;
getFloats(dom, x);
soln[0] = x["temperature"];
resize(1,1);
}
}
/////////////////////////////////////////////////////////////
//
// surf1D

28
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_min_range(0.001),
m_domain(&domain), m_npmax(200)
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_nv = m_domain->nComponents();
m_active.resize(m_nv, true);
@@ -51,6 +51,7 @@ namespace Cantera {
m_keep[0] = 1;
m_keep[n-1] = 1;
//m_did_analysis = false;
if (m_domain->nPoints() <= 1) return 0;
@@ -58,10 +59,13 @@ namespace Cantera {
m_nv = m_domain->nComponents();
// check consistency
if (n != m_domain->nPoints()) throw CanteraError("analyze","inconsistent");
if (n != m_domain->nPoints())
throw CanteraError("analyze","inconsistent");
if (n >= m_npmax) throw CanteraError("analyze","max points");
if (n >= m_npmax) {
return -2; // throw CanteraError("analyze","max points");
}
/**
* find locations where cell size ratio is too large.
@@ -126,12 +130,16 @@ namespace Cantera {
if (r > 1.0) {
m_loc[j] = 1;
m_c[name] = 1;
if (int(m_loc.size()) + n > m_npmax) goto done;
//if (int(m_loc.size()) + n > m_npmax) goto done;
}
if (r >= -1.0) {
if (r >= m_prune) {
m_keep[j] = 1;
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;
}
}
}
@@ -152,17 +160,21 @@ namespace Cantera {
m_c[name] = 1;
m_loc[j] = 1;
m_loc[j+1] = 1;
if (int(m_loc.size()) + n > m_npmax) goto done;
//if (int(m_loc.size()) + n > m_npmax) goto done;
}
if (r >= -1.0) {
if (r >= m_prune) {
m_keep[j+1] = 1;
}
else {
if (m_keep[j+1] == 0) m_keep[j+1] = -1;
}
}
}
}
}
done:
//m_did_analysis = true;
return m_loc.size();
}

8
Cantera/src/oneD/refine.h
View File

@@ -14,8 +14,10 @@ namespace Cantera {
void setCriteria(doublereal ratio = 10.0,
doublereal slope = 0.8,
doublereal curve = 0.8) {
doublereal curve = 0.8,
doublereal prune = 0.1) {
m_ratio = ratio; m_slope = slope; m_curve = curve;
m_prune = prune;
}
void setActive(int comp, bool state = true) { m_active[comp] = state; }
void setMaxPoints(int npmax) { m_npmax = npmax; }
@@ -28,7 +30,7 @@ namespace Cantera {
return m_loc.find(j) != m_loc.end();
}
bool keepPoint(int j) {
return 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);
@@ -38,7 +40,7 @@ namespace Cantera {
map<int, int> m_keep;
map<string, int> m_c;
vector<bool> m_active;
doublereal m_ratio, m_slope, m_curve;
doublereal m_ratio, m_slope, m_curve, m_prune;
doublereal m_min_range;
Domain1D* m_domain;
int m_nv, m_npmax;

1
Cantera/src/xml.cpp
View File

@@ -197,6 +197,7 @@ namespace Cantera {
if (s[s.size()-1] == '/') {
name += "/";
}
// get attributes
while (1) {
iloc = s.find('=');