*** empty log message ***

This commit is contained in:
Dave Goodwin 2005-08-18 14:44:35 +00:00
parent 9cfa1faf5e
commit 83501b5d72
55 changed files with 850 additions and 242 deletions

View file

@ -173,6 +173,23 @@ extern "C" {
return _mix(i)->temperature();
}
doublereal DLL_EXPORT mix_minTemp(int i) {
return _mix(i)->minTemp();
}
doublereal DLL_EXPORT mix_maxTemp(int i) {
return _mix(i)->maxTemp();
}
doublereal DLL_EXPORT mix_charge(int i) {
return _mix(i)->charge();
}
doublereal DLL_EXPORT mix_phaseCharge(int i, int p) {
if (!checkPhase(i,p)) return DERR;
return _mix(i)->phaseCharge(p);
}
int DLL_EXPORT mix_setPressure(int i, double p) {
if (p < 0.0) return -1;
_mix(i)->setPressure(p);
@ -217,4 +234,47 @@ extern "C" {
}
}
int DLL_EXPORT mix_getValidChemPotentials(int i, double bad_mu,
int standard, int lenmu, double* mu) {
bool st = (standard == 1);
try {
if (lenmu < _mix(i)->nSpecies())
throw CanteraError("getChemPotentials","array too small");
_mix(i)->getValidChemPotentials(bad_mu, mu, st);
return 0;
}
catch (CanteraError) {
return -1;
}
}
double DLL_EXPORT mix_enthalpy(int i) {
return _mix(i)->enthalpy();
}
double DLL_EXPORT mix_entropy(int i) {
return _mix(i)->entropy();
}
double DLL_EXPORT mix_gibbs(int i) {
return _mix(i)->gibbs();
}
double DLL_EXPORT mix_cp(int i) {
return _mix(i)->cp();
}
double DLL_EXPORT mix_volume(int i) {
return _mix(i)->volume();
}
int DLL_EXPORT mix_speciesPhaseIndex(int i, int k) {
return _mix(i)->speciesPhaseIndex(k);
}
double DLL_EXPORT mix_moleFraction(int i, int k) {
return _mix(i)->moleFraction(k);
}
}

View file

@ -17,6 +17,10 @@ extern "C" {
int DLL_IMPORT mix_nSpecies(int i);
int DLL_IMPORT mix_setTemperature(int i, double t);
double DLL_IMPORT mix_temperature(int i);
double DLL_IMPORT mix_minTemp(int i);
double DLL_IMPORT mix_maxTemp(int i);
double DLL_IMPORT mix_charge(int i);
double DLL_IMPORT mix_phaseCharge(int i, int p);
int DLL_IMPORT mix_setPressure(int i, double p);
double DLL_IMPORT mix_pressure(int i);
double DLL_IMPORT mix_nAtoms(int i, int k, int m);
@ -30,5 +34,19 @@ extern "C" {
double DLL_IMPORT mix_equilibrate(int i, char* XY,
double err, int maxsteps, int maxiter, int loglevel);
int DLL_IMPORT mix_getChemPotentials(int i, int lenmu, double* mu);
int DLL_IMPORT mix_getValidChemPotentials(int i, double bad_mu,
int standard, int lenmu, double* mu);
double DLL_IMPORT mix_enthalpy(int i);
double DLL_IMPORT mix_entropy(int i);
double DLL_IMPORT mix_gibbs(int i);
double DLL_IMPORT mix_cp(int i);
double DLL_IMPORT mix_volume(int i);
int DLL_IMPORT mix_speciesPhaseIndex(int k);
double DLL_IMPORT mix_moleFraction(int k);
}
#endif

View file

@ -483,7 +483,7 @@ extern "C" {
return 0;
}
int DLL_EXPORT sim1D_setTimeStep(int i, double stepsize, int ns, int* nsteps) {
int DLL_EXPORT sim1D_setTimeStep(int i, double stepsize, int ns, integer* nsteps) {
try {
_sim1D(i)->setTimeStep(stepsize, ns, nsteps);
return 0;

View file

@ -2,6 +2,7 @@
#define CTC_ONEDIM_H
#include "clib_defs.h"
#include "../../src/config.h"
extern "C" {
@ -61,7 +62,7 @@ extern "C" {
int np, double* pos, int nv, double* v);
int DLL_IMPORT sim1D_setFlatProfile(int i, int dom, int comp, double v);
int DLL_IMPORT sim1D_showSolution(int i, char* fname);
int DLL_IMPORT sim1D_setTimeStep(int i, double stepsize, int ns, int* nsteps);
int DLL_IMPORT sim1D_setTimeStep(int i, double stepsize, int ns, integer* nsteps);
int DLL_IMPORT sim1D_getInitialSoln(int i);
int DLL_IMPORT sim1D_solve(int i, int loglevel, int refine_grid);
int DLL_IMPORT sim1D_refine(int i, int loglevel);

View file

@ -27,7 +27,7 @@ FORT_LIBS = @FLIBS@
CXX = @CXX@
# C++ compile flags
CXX_FLAGS = @CXXFLAGS@
CXX_FLAGS = @CXXFLAGS@ @CXX_INCLUDES@
# external libraries
EXT_LIBS = @LOCAL_LIBS@ -lctcxx

View file

@ -84,6 +84,10 @@ class Mixture:
"""Total number of phases defined for the mixture."""
return len(self._phases)
def phase(self, n):
"""Return the object representing the nth phase in the mixture."""
return self._phases[n]
def phaseName(self, n):
"""Name of phase n."""
return self._phases[n].name()
@ -162,16 +166,40 @@ class Mixture:
"""Set the temperature [K]. The temperatures of all phases are
set to this value, holding the pressure fixed."""
return _cantera.mix_setTemperature(self.__mixid, t)
def temperature(self):
"""The temperature [K]."""
return _cantera.mix_temperature(self.__mixid)
def minTemp(self):
"""The minimum temperature for which all species in
multi-species solutions have valid thermo data. Stoichiometric
phases are not considered in determining minTemp. """
return _cantera.mix_minTemp(self._mixid)
def maxTemp(self):
"""The maximum temperature for which all species in
multi-species solutions have valid thermo data. Stoichiometric
phases are not considered in determining maxTemp. """
return _cantera.mix_maxTemp(self._index)
def charge(self):
"""The total charge in Coulombs, summed over all phases."""
return _cantera.mix_charge(self._index)
def phaseCharge(self, p):
"""The charge of phase p (Coulombs)."""
return _cantera.mix_phaseCharge(self._index, p)
def setPressure(self, p):
"""Set the pressure [Pa]. The pressures of all phases are set
to the specified value, holding the temperature fixed."""
return _cantera.mix_setPressure(self.__mixid, p)
def pressure(self):
"""The pressure [Pa]."""
return _cantera.mix_pressure(self.__mixid)
return _cantera.mix_pressure(self.__mixid)
def phaseMoles(self, n = -1):
"""Moles of phase n."""
if n == -1:
@ -182,9 +210,11 @@ class Mixture:
return moles
else:
return _cantera.mix_phaseMoles(self.__mixid, n)
def setPhaseMoles(self, n, moles):
"""Set the number of moles of phase n."""
_cantera.mix_setPhaseMoles(self.__mixid, n, moles)
def setSpeciesMoles(self, moles):
"""Set the moles of the species [kmol]. The moles may be
specified either as a string, or as an array. If an array is

View file

@ -135,6 +135,63 @@ py_mix_setTemperature(PyObject *self, PyObject *args)
}
static PyObject *
py_mix_minTemp(PyObject *self, PyObject *args)
{
double _val;
int i;
if (!PyArg_ParseTuple(args, "i:mix_minTemp", &i))
return NULL;
_val = mix_minTemp(i);
if (int(_val) == -1) return reportCanteraError();
return Py_BuildValue("d",_val);
}
static PyObject *
py_mix_maxTemp(PyObject *self, PyObject *args)
{
double _val;
int i;
if (!PyArg_ParseTuple(args, "i:mix_maxTemp", &i))
return NULL;
_val = mix_maxTemp(i);
if (int(_val) == -1) return reportCanteraError();
return Py_BuildValue("d",_val);
}
static PyObject *
py_mix_charge(PyObject *self, PyObject *args)
{
double _val;
int i;
if (!PyArg_ParseTuple(args, "i:mix_charge", &i))
return NULL;
_val = mix_charge(i);
if (int(_val) == -1) return reportCanteraError();
return Py_BuildValue("d",_val);
}
static PyObject *
py_mix_phaseCharge(PyObject *self, PyObject *args)
{
double _val;
int i;
int p;
if (!PyArg_ParseTuple(args, "ii:mix_phaseCharge", &i, &p))
return NULL;
_val = mix_phaseCharge(i,p);
if (int(_val) == -1) return reportCanteraError();
return Py_BuildValue("d",_val);
}
static PyObject *
py_mix_temperature(PyObject *self, PyObject *args)
{

View file

@ -713,7 +713,7 @@ py_sim1D_setTimeStep(PyObject *self, PyObject *args)
PyArrayObject* nsteps_array = (PyArrayObject*)
PyArray_ContiguousFromObject(nsteps, PyArray_INT, 1, 1);
int* nsteps_data = (int*)nsteps_array->data;
integer* nsteps_data = (integer*)nsteps_array->data;
int nsteps_len = nsteps_array->dimensions[0];
_val = sim1D_setTimeStep(i,stepsize,nsteps_len,nsteps_data);

View file

@ -265,6 +265,10 @@ static PyMethodDef ct_methods[] = {
{"mix_nAtoms", py_mix_nAtoms, METH_VARARGS},
{"mix_setTemperature", py_mix_setTemperature, METH_VARARGS},
{"mix_temperature", py_mix_temperature, METH_VARARGS},
{"mix_minTemp", py_mix_minTemp, METH_VARARGS},
{"mix_maxTemp", py_mix_maxTemp, METH_VARARGS},
{"mix_charge", py_mix_charge, METH_VARARGS},
{"mix_phaseCharge", py_mix_phaseCharge, METH_VARARGS},
{"mix_setPressure", py_mix_setPressure, METH_VARARGS},
{"mix_pressure", py_mix_pressure, METH_VARARGS},
{"mix_phaseMoles", py_mix_phaseMoles, METH_VARARGS},

View file

@ -64,7 +64,7 @@ namespace Cantera {
void ConstDensityThermo::getActivityCoefficients(doublereal* ac) const {
for (int k = 0; k < m_kk; k++) {
ac[k] = 1.0;
ac[k] = 1.0;
}
}
@ -132,7 +132,7 @@ namespace Cantera {
}
void ConstDensityThermo::setParametersFromXML(const XML_Node& eosdata) {
eosdata.require("model","Incompressible");
eosdata._require("model","Incompressible");
doublereal rho = getFloat(eosdata, "density", "-");
setDensity(rho);
}

View file

@ -108,7 +108,7 @@ namespace Cantera {
if (info != 0)
throw CanteraError("DenseMatrix::leaseSquares",
"DGELSS returned INFO = "+int2str(info));
return 0;
return 0;
}
#endif

View file

@ -21,6 +21,8 @@
#include "ThirdBodyMgr.h"
#include "RateCoeffMgr.h"
//#include "../user/grirxnstoich.h"
#include <iostream>
using namespace std;
@ -44,11 +46,11 @@ namespace Cantera {
if (thermo != 0) addPhase(*thermo);
m_kdata = new GasKineticsData;
m_kdata->m_temp = 0.0;
// m_rxnstoich = new ReactionStoichMgr;
m_rxnstoich = new ReactionStoichMgr;
}
GasKinetics::
~GasKinetics() {delete m_kdata;}
~GasKinetics() {delete m_kdata; delete m_rxnstoich;}
/**
* Update temperature-dependent portions of reaction rates and
@ -103,7 +105,7 @@ namespace Cantera {
fill(m_rkc.begin(), m_rkc.end(), 0.0);
// compute Delta G^0 for all reversible reactions
m_rxnstoich.getRevReactionDelta(m_ii, m_grt.begin(), m_rkc.begin());
m_rxnstoich->getRevReactionDelta(m_ii, m_grt.begin(), m_rkc.begin());
doublereal logStandConc = m_kdata->m_logStandConc;
doublereal rrt = 1.0/(GasConstant * thermo().temperature());
@ -130,7 +132,7 @@ namespace Cantera {
fill(rkc.begin(), rkc.end(), 0.0);
// compute Delta G^0 for all reactions
m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), rkc.begin());
m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), rkc.begin());
doublereal logStandConc = m_kdata->m_logStandConc;
doublereal rrt = 1.0/(GasConstant * thermo().temperature());
@ -160,7 +162,7 @@ namespace Cantera {
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaG);
m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaG);
}
/**
@ -184,7 +186,7 @@ namespace Cantera {
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaH);
m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaH);
}
/************************************************************************
@ -208,7 +210,7 @@ namespace Cantera {
* Use the stoichiometric manager to find deltaS for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaS);
m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaS);
}
/**
@ -234,7 +236,7 @@ namespace Cantera {
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaG);
m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaG);
}
/**
@ -264,7 +266,7 @@ namespace Cantera {
* Use the stoichiometric manager to find deltaG for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaH);
m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaH);
}
/*********************************************************************
@ -293,7 +295,7 @@ namespace Cantera {
* Use the stoichiometric manager to find deltaS for each
* reaction.
*/
m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaS);
m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaS);
}
void GasKinetics::processFalloffReactions() {
@ -356,12 +358,12 @@ namespace Cantera {
multiply_each(ropr.begin(), ropr.end(), m_rkc.begin());
// multiply ropf by concentration products
m_rxnstoich.multiplyReactants(m_conc.begin(), ropf.begin());
m_rxnstoich->multiplyReactants(m_conc.begin(), ropf.begin());
//m_reactantStoich.multiply(m_conc.begin(), ropf.begin());
// for reversible reactions, multiply ropr by concentration
// products
m_rxnstoich.multiplyRevProducts(m_conc.begin(), ropr.begin());
m_rxnstoich->multiplyRevProducts(m_conc.begin(), ropr.begin());
//m_revProductStoich.multiply(m_conc.begin(), ropr.begin());
for (int j = 0; j != m_ii; ++j) {
@ -578,7 +580,7 @@ namespace Cantera {
m_kdata->m_rkcn.push_back(0.0);
m_rxnstoich.add(reactionNumber(), r);
m_rxnstoich->add(reactionNumber(), r);
if (r.reversible) {
m_dn.push_back(pk.size() - rk.size());
@ -632,6 +634,7 @@ namespace Cantera {
// m_pstoich[i][m_products[i][j]]++;
// }
// }
//m_rxnstoich->write("c.cpp");
m_finalized = true;
}
}

View file

@ -58,7 +58,7 @@ namespace Cantera {
bool m_ROP_ok;
doublereal m_temp;
vector_fp m_rfn;
vector_fp m_rfn;
vector_fp falloff_work;
vector_fp concm_3b_values;
vector_fp concm_falloff_values;
@ -218,7 +218,7 @@ namespace Cantera {
#ifdef HWMECH
get_wdot(m_kdata->m_ropnet.begin(), net);
#else
m_rxnstoich.getNetProductionRates(m_kk, m_kdata->m_ropnet.begin(), net);
m_rxnstoich->getNetProductionRates(m_kk, m_kdata->m_ropnet.begin(), net);
//fill(net, net + m_kk, 0.0);
//m_revProductStoich.incrementSpecies(
// m_kdata->m_ropnet.begin(), net);
@ -238,7 +238,7 @@ namespace Cantera {
*/
virtual void getCreationRates(doublereal* cdot) {
updateROP();
m_rxnstoich.getCreationRates(m_kk, m_kdata->m_ropf.begin(),
m_rxnstoich->getCreationRates(m_kk, m_kdata->m_ropf.begin(),
m_kdata->m_ropr.begin(), cdot);
//fill(cdot, cdot + m_kk, 0.0);
//m_revProductStoich.incrementSpecies(
@ -258,7 +258,7 @@ namespace Cantera {
*/
virtual void getDestructionRates(doublereal* ddot) {
updateROP();
m_rxnstoich.getDestructionRates(m_kk, m_kdata->m_ropf.begin(),
m_rxnstoich->getDestructionRates(m_kk, m_kdata->m_ropf.begin(),
m_kdata->m_ropr.begin(), ddot);
// fill(ddot, ddot + m_kk, 0.0);
//m_revProductStoich.incrementSpecies(
@ -380,7 +380,7 @@ namespace Cantera {
//StoichManagerN m_revProductStoich;
//StoichManagerN m_irrevProductStoich;
ReactionStoichMgr m_rxnstoich;
ReactionStoichMgr* m_rxnstoich;
vector<int> m_fwdOrder;

View file

@ -116,7 +116,7 @@ namespace Cantera {
}
void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata) {
eosdata.require("model","LatticeSolid");
eosdata._require("model","LatticeSolid");
XML_Node& la = eosdata.child("LatticeArray");
vector<XML_Node*> lattices;
la.getChildren("Lattice",lattices);

View file

@ -160,7 +160,7 @@ flow1D:
cd oneD; @MAKE@
CXX_LIBS = @LIBS@
CXX_INCLUDES = -I.
CXX_INCLUDES = @CXX_INCLUDES@ -I.
CANTERA_LIB = @buildlib@/libcantera.a
DEPENDS = $(ALL_OBJ:.o=.d)
@ -172,7 +172,7 @@ ALL_H = $(BASE_H) $(THERMO_H) $(KINETICS_H) $(HETEROKIN_H) \
g++ -MM $(CXX_INCLUDES) $*.cpp > $*.d
.cpp.o:
@CXX@ -c $< $(CXX_FLAGS)
@CXX@ -c $< $(CXX_INCLUDES) $(CXX_FLAGS)
lib: $(OBJ_LIB)
$(RM) $(CANTERA_LIB)

View file

@ -75,7 +75,7 @@ namespace Cantera {
}
virtual void setParametersFromXML(const XML_Node& eosdata) {
eosdata.require("model","Metal");
eosdata._require("model","Metal");
doublereal rho = getFloat(eosdata, "density", "-");
setDensity(rho);
}

View file

@ -441,7 +441,11 @@ namespace Cantera {
for (n = 0; n < maxiter; n++) {
// if 'strt' is false, the current composition will be used as
// the starting estimate; otherwise it will be estimated
// the starting estimate; otherwise it will be estimated
// if (e) {
// cout << "e should be zero, but it is not!" << endl;
// delete e;
// }
e = new MultiPhaseEquil(this, strt);
// start with a loose error tolerance, but tighten it as we get
// close to the final temperature
@ -521,9 +525,9 @@ namespace Cantera {
}
endLogGroup();
}
delete e;
e = 0;
}
delete e;
e = 0;
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
throw CanteraError("MultiPhase::equilibrate",
@ -541,6 +545,7 @@ namespace Cantera {
addLogEntry("max T",fp2str(Thigh));
for (n = 0; n < maxiter; n++) {
if (e) delete e;
e = new MultiPhaseEquil(this, strt);
ferr = 0.1;
if (fabs(dt) < 1.0) ferr = err;
@ -596,9 +601,9 @@ namespace Cantera {
}
endLogGroup();
}
delete e;
e = 0;
}
delete e;
e = 0;
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
throw CanteraError("MultiPhase::equilibrate",

View file

@ -102,13 +102,13 @@ namespace Cantera {
/// valid thermo data. Stoichiometric phases are not
/// considered, since they may have thermo data only valid for
/// conditions for which they are stable.
doublereal minTemp();
doublereal minTemp() { return m_Tmin; }
/// Maximum temperature for which all solution phases have
/// valid thermo data. Stoichiometric phases are not
/// considered, since they may have thermo data only valid for
/// conditions for which they are stable.
doublereal maxTemp();
doublereal maxTemp() { return m_Tmax; }
/// Total charge (Coulombs).
doublereal charge();

View file

@ -228,6 +228,8 @@ namespace Cantera {
index_t m, n, ik, j;
double not_mu = 1.0e12;
beginLogGroup("MultiPhaseEquil::setInitialMoles");
m_mix->getValidChemPotentials(not_mu, m_mu.begin(), true);
doublereal dg_rt;
@ -241,7 +243,7 @@ namespace Cantera {
// choose a set of components based on the current
// composition
computeN();
addLogEntry("iteration",iter);
redo = false;
iter++;
if (iter > 4) break;
@ -265,7 +267,10 @@ namespace Cantera {
delta_xi = fabs(moles(ik)/nu);
// if a component has nearly zero moles, redo
// with a new set of components
if (delta_xi < SmallNumber && ik < m_nel) redo = true;
if (!redo && delta_xi < 1.0e-10 && ik < m_nel) {
addLogEntry("component too small",speciesName(ik));
redo = true;
}
if (delta_xi < dxi_min) dxi_min = delta_xi;
}
}
@ -277,7 +282,10 @@ namespace Cantera {
// set the moles of the phase objects to match
updateMixMoles();
}
for (ik = 0; ik < m_nsp; ik++)
if (moles(ik) != 0.0) addLogEntry(speciesName(ik), moles(ik));
endLogGroup("MultiPhaseEquil::setInitialMoles");
return 0;
}
@ -446,21 +454,21 @@ namespace Cantera {
k = m_species[ik];
addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik]));
}
endLogGroup();
endLogGroup("components");
beginLogGroup("non-components");
for (m = m_nel; m < m_nsp; m++) {
ik = m_order[m];
k = m_species[ik];
addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik]));
}
endLogGroup();
endLogGroup("non-components");
addLogEntry("Error",fp2str(error()));
beginLogGroup("Delta G / RT");
for (k = 0; k < m_nsp - m_nel; k++) {
addLogEntry(reactionString(k), fp2str(m_deltaG_RT[k]));
}
endLogGroup();
endLogGroup();
endLogGroup("Delta G / RT");
endLogGroup("info");
}
/// Return a string specifying the jth reaction.
@ -582,7 +590,6 @@ namespace Cantera {
// now take a step with this scaled omega
addLogEntry("Stepping by ", fp2str(omegamax));
step(omegamax, m_work);
// compute the gradient of G at this new position in the
// current direction. If it is positive, then we have overshot
// the minimum. In this case, interpolate back.
@ -600,6 +607,7 @@ namespace Cantera {
addLogEntry("Stepped over minimum. Take smaller step ", fp2str(omega));
step(omega, m_work);
}
printInfo();
endLogGroup("MultiPhaseEquil::stepComposition");
return omega;
}

View file

@ -44,7 +44,7 @@ namespace Cantera {
void PureFluidPhase::
setParametersFromXML(const XML_Node& eosdata) {
eosdata.require("model","PureFluid");
eosdata._require("model","PureFluid");
m_subflag = atoi(eosdata["fluid_type"].c_str());
if (m_subflag < 0)
throw CanteraError("PureFluidPhase::setParametersFromXML",

View file

@ -32,7 +32,9 @@ namespace Cantera {
m_reactants = new StoichManagerN;
m_revproducts = new StoichManagerN;
m_irrevproducts = new StoichManagerN;
//m_global = new StoichManagerN;
#ifdef INCL_STOICH_WRITER
m_rwriter = new StoichWriter;
#endif
m_dummy.resize(10,1.0);
}
@ -42,6 +44,9 @@ namespace Cantera {
delete m_revproducts;
delete m_irrevproducts;
// delete m_global;
#ifdef INCL_STOICH_WRITER
delete m_rwriter;
#endif
}
@ -78,12 +83,17 @@ namespace Cantera {
// or specified reaction orders, then add it in a ma
if (isfrac || r.global || rk.size() > 3) {
m_reactants->add(rxn, r.reactants, r.order, r.rstoich);
#ifdef INCL_STOICH_WRITER
if (m_rwriter) m_rwriter->add(rxn, r.reactants, r.order, r.rstoich);
#endif
}
else {
m_reactants->add( rxn, rk);
#ifdef INCL_STOICH_WRITER
if (m_rwriter) m_rwriter->add(rxn, rk);
#endif
}
vector_int pk;
isfrac = false;
int np = r.products.size();
@ -177,4 +187,93 @@ namespace Cantera {
multiplyRevProducts(const doublereal* c, doublereal* r) {
m_revproducts->multiply(c, r);
}
void ReactionStoichMgr::
write(string filename) {
ofstream f(filename.c_str());
f << "namespace mech {" << endl;
writeCreationRates(f);
writeDestructionRates(f);
writeNetProductionRates(f);
writeMultiplyReactants(f);
writeMultiplyRevProducts(f);
f << "} // namespace mech" << endl;
f.close();
}
void ReactionStoichMgr::
writeCreationRates(ostream& f) {
f << " void getCreationRates(const doublereal* rf, const doublereal* rb," << endl;
f << " doublereal* c) {" << endl;
map<int, string> out;
m_revproducts->writeIncrementSpecies("rf",out);
m_irrevproducts->writeIncrementSpecies("rf",out);
m_reactants->writeIncrementSpecies("rb",out);
map<int, string>::iterator b;
for (b = out.begin(); b != out.end(); ++b) {
string rhs = wrapString(b->second);
rhs[1] = '=';
f << " c[" << b->first << "] " << rhs << ";" << endl;
}
f << " }" << endl << endl << endl;
}
void ReactionStoichMgr::
writeDestructionRates(ostream& f) {
f << " void getDestructionRates(const doublereal* rf, const doublereal* rb," << endl;
f << " doublereal* d) {" << endl;
map<int, string> out;
m_revproducts->writeIncrementSpecies("rb",out);
m_reactants->writeIncrementSpecies("rf",out);
map<int, string>::iterator b;
for (b = out.begin(); b != out.end(); ++b) {
string rhs = wrapString(b->second);
rhs[1] = '=';
f << " d[" << b->first << "] " << rhs << ";" << endl;
}
f << " }" << endl << endl << endl;
}
void ReactionStoichMgr::
writeNetProductionRates(ostream& f) {
f << " void getNetProductionRates(const doublereal* r, doublereal* w) {" << endl;
map<int, string> out;
m_revproducts->writeIncrementSpecies("r",out);
m_irrevproducts->writeIncrementSpecies("r",out);
m_reactants->writeDecrementSpecies("r",out);
map<int, string>::iterator b;
for (b = out.begin(); b != out.end(); ++b) {
string rhs = wrapString(b->second);
rhs[1] = '=';
f << " w[" << b->first << "] " << rhs << ";" << endl;
}
f << " }" << endl << endl << endl;
}
void ReactionStoichMgr::
writeMultiplyReactants(ostream& f) {
f << " void multiplyReactants(const doublereal* c, doublereal* r) {" << endl;
map<int, string> out;
m_reactants->writeMultiply("c",out);
map<int, string>::iterator b;
for (b = out.begin(); b != out.end(); ++b) {
string rhs = b->second;
f << " r[" << b->first << "] *= " << rhs << ";" << endl;
}
f << " }" << endl << endl << endl;
}
void ReactionStoichMgr::
writeMultiplyRevProducts(ostream& f) {
f << " void multiplyRevProducts(const doublereal* c, doublereal* r) {" << endl;
map<int, string> out;
m_revproducts->writeMultiply("c",out);
map<int, string>::iterator b;
for (b = out.begin(); b != out.end(); ++b) {
string rhs = b->second;
f << " r[" << b->first << "] *= " << rhs << ";" << endl;
}
f << " }" << endl << endl << endl;
}
}

View file

@ -96,7 +96,7 @@ namespace Cantera {
* @param products vector of integer product indices
* @param reversible true if the reaction is reversible, false otherwise
*/
void add(int rxn, const vector_int& reactants, const vector_int& products,
virtual void add(int rxn, const vector_int& reactants, const vector_int& products,
bool reversible);
/**
@ -116,7 +116,7 @@ namespace Cantera {
// bool reversible, const vector_fp& fwdOrder);
void add(int rxn, const ReactionData& r);
virtual void add(int rxn, const ReactionData& r);
/**
* Species creation rates.
@ -127,7 +127,7 @@ namespace Cantera {
* C = N_p Q_f + N_r Q_r.
* \f]
*/
void getCreationRates(int nSpecies,
virtual void getCreationRates(int nSpecies,
const doublereal* fwdRatesOfProgress,
const doublereal* revRatesOfProgress,
doublereal* creationRates);
@ -144,7 +144,7 @@ namespace Cantera {
* Note that the stoichiometric coefficient matrices are very sparse, integer
* matrices.
*/
void getDestructionRates(int nSpecies,
virtual void getDestructionRates(int nSpecies,
const doublereal* fwdRatesOfProgress,
const doublereal* revRatesOfProgress,
doublereal* destructionRates);
@ -164,7 +164,7 @@ namespace Cantera {
* W = (N_r - N_p) Q_{\rm net},
* \f]
*/
void getNetProductionRates(int nsp, const doublereal* ropnet, doublereal* w);
virtual void getNetProductionRates(int nsp, const doublereal* ropnet, doublereal* w);
@ -176,7 +176,7 @@ namespace Cantera {
* and array 'dg' must have a length as great as the total
* number of reactions.
*/
void getReactionDelta(int nReactions,
virtual void getReactionDelta(int nReactions,
const doublereal* g,
doublereal* dg);
@ -193,7 +193,7 @@ namespace Cantera {
* calculating reveerse rate coefficients from thermochemistry
* for reversible reactions.
*/
void getRevReactionDelta(int nr, const doublereal* g, doublereal* dg);
virtual void getRevReactionDelta(int nr, const doublereal* g, doublereal* dg);
/**
@ -204,7 +204,7 @@ namespace Cantera {
* \f]
* Here \f$ o_{k,i} \f$ is the reaction order of species k in reaction i.
*/
void multiplyReactants(const doublereal* C, doublereal* R);
virtual void multiplyReactants(const doublereal* C, doublereal* R);
/**
@ -216,15 +216,25 @@ namespace Cantera {
* Here \f$ \nu^{(p)}_{k,i} \f$ is the product stoichiometric coefficient
* of species k in reaction i.
*/
void multiplyRevProducts(const doublereal* c, doublereal* r);
virtual void multiplyRevProducts(const doublereal* c, doublereal* r);
virtual void write(string filename);
protected:
void writeCreationRates(ostream& f);
void writeDestructionRates(ostream& f);
void writeNetProductionRates(ostream& f);
void writeMultiplyReactants(ostream& f);
void writeMultiplyRevProducts(ostream& f);
StoichManagerN* m_reactants;
StoichManagerN* m_revproducts;
StoichManagerN* m_irrevproducts;
vector_fp m_dummy;
#ifdef INCL_STOICH_WRITER
StoichWriter* m_rwriter;
#endif
};
}

View file

@ -19,7 +19,9 @@
#include "ctexceptions.h"
#include "stringUtils.h"
#include "State.h"
#ifdef DARWIN
#include <Accelerate.h>
#endif
namespace Cantera {
State::State() : m_kk(0), m_temp(0.0), m_dens(0.001), m_mmw(0.0) {}
@ -71,10 +73,14 @@ namespace Cantera {
void State::setMassFractions(const doublereal* y) {
doublereal norm = 0.0, sum = 0.0;
int k;
cblas_dcopy(m_kk, y, 1, m_y.begin(), 1);
for (k = 0; k != m_kk; ++k) {
norm += y[k];
//m_y[k] = y[k];
}
scale(y, y + m_kk, m_y.begin(), 1.0/norm);
//scale(y, y + m_kk, m_y.begin(), 1.0/norm);
scale(m_kk, 1.0/norm, m_y.begin());
for (k = 0; k != m_kk; ++k) {
m_ym[k] = m_y[k] * m_rmolwts[k];
sum += m_ym[k];
@ -118,6 +124,7 @@ namespace Cantera {
}
void State::getConcentrations(doublereal* c) const {
//ct_dscal(m_kk, m_dens, m_ym.begin(), 1);
scale(m_ym.begin(), m_ym.end(), c, m_dens);
}
@ -130,6 +137,7 @@ namespace Cantera {
}
void State::getMoleFractions(doublereal* x) const {
//ct_dscal(m_kk, m_mmw, m_ym.begin(), 1);
scale(m_ym.begin(), m_ym.end(), x, m_mmw);
}
@ -161,27 +169,3 @@ namespace Cantera {
}
}

View file

@ -22,8 +22,55 @@ namespace Cantera {
* Note: these classes are designed for internal use in class
* ReactionStoichManager.
*
* The classes defined here implement simple operations that are
* used by class ReactionStoichManager to compute things like
* rates of progress, species production rates, etc. In general, a
* reaction mechanism may involve many species and many reactions,
* but any given reaction typically only involves a few species as
* reactants, and a few as products. Therefore, the matrix of
* stoichiometric coefficients is very sparse. Not only is it
* sparse, but the non-zero matrix elements often have the value
* 1, and in many cases no more than three coefficients are
* non-zero for the reactants and/or the products.
*
* For the present purposes, we will consider each direction of a
* reversible reaction to be a separate reaction. We often need to
* compute quantities that can formally be written as a matrix
* product of a stoichiometric coefficient matrix and a vector of
* reaction rates. For example, the species creation rates are
* given by
* \f[
* \dot C_k = \sum_k \nu^{(p)}_{k,i} R_i
* \f]
* where \f$ \nu^{(p)_{k,i}}$ is the product-side stoichiometric
* coefficient of species \a k in reaction \a i.
* This could be done be straightforward matrix multiplication, but would be inefficient, since most of the matrix elements of \f$ \nu^{(p)}_{k,i} \f$ are zero. We could do better by using sparse-matrix algorithms to compute this product.
If the reactions are general ones, with non-integral stoichiometric
coefficients, this is about as good as we can do. But we are
particularly concerned here with the performance for very large
reaction mechanisms, which are usually composed of elementary
reactions, which have integral stoichiometric
coefficients. Furthermore, very few elementary reactions involve more
than 3 product or reactant molecules. This means that instead of
But we can do even better if we take account of the special structure
of this matrix for elementary reactions.
involve three or fewer product molecules (or reactant molecules).
* To take advantage of this structure, reactions are divided int
These classes are
* designed to take advantage of this sparse structure when
* computing quantities that can be written as matrix multiplies
They are designed to explicitly unroll loops over species or reactions for
* Operations on reactions that require knowing the reaction
* stoichiometry. This module consists of class StoichManager, and
* stoichiometry.
* This module consists of class StoichManager, and
* classes C1, C2, and C3. Classes C1, C2, and C3 handle operations
* involving one, two, or three species, respectively, in a
* reaction. Instances are instantiated with a reaction number, and n
@ -101,6 +148,9 @@ namespace Cantera {
return 0.0;
}
inline static string fmt(string r, int n) { return r + "[" + int2str(n) + "]"; }
/**
* Handles one species in a reaction.
* @ingroup Stoichiometry
@ -139,6 +189,28 @@ namespace Cantera {
R[m_rxn] -= S[m_ic0];
}
int rxnNumber() const { return m_rxn; }
int speciesIndex(int n) const { return m_ic0; }
int nSpecies() { return 1;}
void writeMultiply(string r, map<int, string>& out) {
out[m_rxn] = fmt(r, m_ic0);
}
void writeIncrementReaction(string r, map<int, string>& out) {
out[m_rxn] += " + "+fmt(r, m_ic0);
}
void writeDecrementReaction(string r, map<int, string>& out) {
out[m_rxn] += " - "+fmt(r, m_ic0);
}
void writeIncrementSpecies(string r, map<int, string>& out) {
out[m_ic0] += " + "+fmt(r, m_rxn);
}
void writeDecrementSpecies(string r, map<int, string>& out) {
out[m_ic0] += " - "+fmt(r, m_rxn);
}
private:
int m_rxn, m_ic0;
};
@ -183,33 +255,30 @@ namespace Cantera {
R[m_rxn] -= (S[m_ic0] + S[m_ic1]);
}
// void multiply(const doublereal* input, doublereal* output) const {
// output[m_rxn] *= input[m_ic0] * input[m_ic1];
// }
int rxnNumber() const { return m_rxn; }
int speciesIndex(int n) const { return (n == 0 ? m_ic0 : m_ic1); }
int nSpecies() { return 2;}
// void incrementSpecies(const doublereal* input,
// doublereal* output) const {
// doublereal x = input[m_rxn];
// output[m_ic0] += x;
// output[m_ic1] += x;
// }
void writeMultiply(string r, map<int, string>& out) {
out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1);
}
void writeIncrementReaction(string r, map<int, string>& out) {
out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1);
}
void writeDecrementReaction(string r, map<int, string>& out) {
out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1);
}
// void decrementSpecies(const doublereal* input,
// doublereal* output) const {
// doublereal x = input[m_rxn];
// output[m_ic0] -= x;
// output[m_ic1] -= x;
// }
// void incrementReaction(const doublereal* input,
// doublereal* output) const {
// *(output + m_rxn) += *(input + m_ic0) + *(input + m_ic1);
// }
// void decrementReaction(const doublereal* input,
// doublereal* output) const {
// *(output + m_rxn) -= (*(input + m_ic0) + *(input + m_ic1));
// }
void writeIncrementSpecies(string r, map<int, string>& out) {
string s = " + "+fmt(r, m_rxn);
out[m_ic0] += s;
out[m_ic1] += s;
}
void writeDecrementSpecies(string r, map<int, string>& out) {
string s = " - "+fmt(r, m_rxn);
out[m_ic0] += s;
out[m_ic1] += s;
}
private:
@ -267,35 +336,31 @@ namespace Cantera {
R[m_rxn] -= (S[m_ic0] + S[m_ic1] + S[m_ic2]);
}
// void multiply(const doublereal* input, doublereal* output) const {
// *(output + m_rxn) *= (*(input + m_ic0)) * (*(input + m_ic1))
// * (*(input + m_ic2));
// }
// void incrementSpecies(const doublereal* input,
// doublereal* output) const {
// doublereal x = *(input + m_rxn);
// *(output + m_ic0) += x;
// *(output + m_ic1) += x;
// *(output + m_ic2) += x;
// }
// void decrementSpecies(const doublereal* input,
// doublereal* output) const {
// doublereal x = *(input + m_rxn);
// *(output + m_ic0) -= x;
// *(output + m_ic1) -= x;
// *(output + m_ic2) -= x;
// }
// void incrementReaction(const doublereal* input,
// doublereal* output) const {
// *(output + m_rxn) += *(input + m_ic0) + *(input + m_ic1)
// + *(input + m_ic2);
// }
// void decrementReaction(const doublereal* input,
// doublereal* output) const {
// *(output + m_rxn) -= (*(input + m_ic0) + *(input + m_ic1)
// + *(input + m_ic2));
// }
int rxnNumber() const { return m_rxn; }
int speciesIndex(int n) const { return (n == 0 ? m_ic0 : (n == 1 ? m_ic1 : m_ic2)); }
int nSpecies() { return 3;}
void writeMultiply(string r, map<int, string>& out) {
out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1) + " * " + fmt(r, m_ic2);
}
void writeIncrementReaction(string r, map<int, string>& out) {
out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1)+" + "+fmt(r, m_ic2);
}
void writeDecrementReaction(string r, map<int, string>& out) {
out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1)+" - "+fmt(r, m_ic2);
}
void writeIncrementSpecies(string r, map<int, string>& out) {
string s = " + "+fmt(r, m_rxn);
out[m_ic0] += s;
out[m_ic1] += s;
out[m_ic2] += s;
}
void writeDecrementSpecies(string r, map<int, string>& out) {
string s = " - "+fmt(r, m_rxn);
out[m_ic0] += s;
out[m_ic1] += s;
out[m_ic2] += s;
}
private:
int m_rxn, m_ic0, m_ic1, m_ic2;
};
@ -331,10 +396,9 @@ namespace Cantera {
return m_rxn;
}
//void power(const doublereal* input, doublereal* output) const {
// for (int n = 0; n < m_n; n++) output[m_rxn]
// *= ppow(input[m_ic[n]],m_order[n]);
//}
doublereal order(int n) const {return m_order[n];}
doublereal stoich(int n) const {return m_stoich[n];}
int speciesIndex(int n) const {return m_ic[n];}
void multiply(const doublereal* input, doublereal* output) const {
for (int n = 0; n < m_n; n++) output[m_rxn] *=
@ -365,6 +429,46 @@ namespace Cantera {
-= m_stoich[n]*input[m_ic[n]];
}
void writeMultiply(string r, map<int, string>& out) {
int n;
out[m_rxn] = "";
for (n = 0; n < m_n; n++) {
if (m_order[n] == 1.0)
out[m_rxn] += fmt(r, m_ic[n]);
else
out[m_rxn] += "pow("+fmt(r, m_ic[n])+","+fp2str(m_order[n])+")";
if (n < m_n-1)
out[m_rxn] += " * ";
}
}
void writeIncrementReaction(string r, map<int, string>& out) {
int n;
for (n = 0; n < m_n; n++) {
out[m_rxn] += " + "+fp2str(m_stoich[n]) + "*" + fmt(r, m_ic[n]);
}
}
void writeDecrementReaction(string r, map<int, string>& out) {
int n;
for (n = 0; n < m_n; n++) {
out[m_rxn] += " - "+fp2str(m_stoich[n]) + "*" + fmt(r, m_ic[n]);
}
}
void writeIncrementSpecies(string r, map<int, string>& out) {
string s = fmt(r, m_rxn);
int n;
for (n = 0; n < m_n; n++) {
out[m_ic[n]] += " + "+fp2str(m_stoich[n]) + "*" + s;
}
}
void writeDecrementSpecies(string r, map<int, string>& out) {
string s = fmt(r, m_rxn);
int n;
for (n = 0; n < m_n; n++) {
out[m_ic[n]] += " - "+fp2str(m_stoich[n]) + "*" + s;
}
}
private:
int m_n, m_rxn;
vector_int m_ic;
@ -409,6 +513,36 @@ namespace Cantera {
}
template<class _InputIter>
inline static void _writeIncrementSpecies(_InputIter __begin, _InputIter __end, string r,
map<int, string>& out) {
for (; __begin != __end; ++__begin) __begin->writeIncrementSpecies(r, out);
}
template<class _InputIter>
inline static void _writeDecrementSpecies(_InputIter __begin, _InputIter __end, string r,
map<int, string>& out) {
for (; __begin != __end; ++__begin) __begin->writeDecrementSpecies(r, out);
}
template<class _InputIter>
inline static void _writeIncrementReaction(_InputIter __begin, _InputIter __end, string r,
map<int, string>& out) {
for (; __begin != __end; ++__begin) __begin->writeIncrementReaction(r, out);
}
template<class _InputIter>
inline static void _writeDecrementReaction(_InputIter __begin, _InputIter __end, string r,
map<int, string>& out) {
for (; __begin != __end; ++__begin) __begin->writeDecrementReaction(r, out);
}
template<class _InputIter>
inline static void _writeMultiply(_InputIter __begin, _InputIter __end, string r,
map<int, string>& out) {
for (; __begin != __end; ++__begin) __begin->writeMultiply(r, out);
}
/*
* This class handles operations involving the stoichiometric
* coefficients on one side of a reaction (reactant or product) for
@ -572,6 +706,42 @@ namespace Cantera {
_decrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output);
}
void writeIncrementSpecies(string r, map<int, string>& out) {
_writeIncrementSpecies(m_c1_list.begin(), m_c1_list.end(), r, out);
_writeIncrementSpecies(m_c2_list.begin(), m_c2_list.end(), r, out);
_writeIncrementSpecies(m_c3_list.begin(), m_c3_list.end(), r, out);
_writeIncrementSpecies(m_cn_list.begin(), m_cn_list.end(), r, out);
}
void writeDecrementSpecies(string r, map<int, string>& out) {
_writeDecrementSpecies(m_c1_list.begin(), m_c1_list.end(), r, out);
_writeDecrementSpecies(m_c2_list.begin(), m_c2_list.end(), r, out);
_writeDecrementSpecies(m_c3_list.begin(), m_c3_list.end(), r, out);
_writeDecrementSpecies(m_cn_list.begin(), m_cn_list.end(), r, out);
}
void writeIncrementReaction(string r, map<int, string>& out) {
_writeIncrementReaction(m_c1_list.begin(), m_c1_list.end(), r, out);
_writeIncrementReaction(m_c2_list.begin(), m_c2_list.end(), r, out);
_writeIncrementReaction(m_c3_list.begin(), m_c3_list.end(), r, out);
_writeIncrementReaction(m_cn_list.begin(), m_cn_list.end(), r, out);
}
void writeDecrementReaction(string r, map<int, string>& out) {
_writeDecrementReaction(m_c1_list.begin(), m_c1_list.end(), r, out);
_writeDecrementReaction(m_c2_list.begin(), m_c2_list.end(), r, out);
_writeDecrementReaction(m_c3_list.begin(), m_c3_list.end(), r, out);
_writeDecrementReaction(m_cn_list.begin(), m_cn_list.end(), r, out);
}
void writeMultiply(string r, map<int, string>& out) {
_writeMultiply(m_c1_list.begin(), m_c1_list.end(), r, out);
_writeMultiply(m_c2_list.begin(), m_c2_list.end(), r, out);
_writeMultiply(m_c3_list.begin(), m_c3_list.end(), r, out);
_writeMultiply(m_cn_list.begin(), m_cn_list.end(), r, out);
}
private:
vector<C1> m_c1_list;
@ -590,7 +760,7 @@ namespace Cantera {
map<int, int> m_loc;
};
#undef INCL_STOICH_WRITER
#ifdef INCL_STOICH_WRITER
class StoichWriter {
@ -610,15 +780,33 @@ namespace Cantera {
}
}
void writeIncSpec(ostream& s, int nsp) {
int k;
for (k = 0; k < nsp; k++) {
s << "out[" << k << "] = " << m_is[k] << ";" << endl;
void add(int rxn, const vector_int& k, const vector_fp& order,
const vector_fp& stoich) {
int n, nn = k.size();
string s;
for (n = 0; n < nn; n++) {
if (order[n] == 1.0)
m_mult[rxn] += "*c[" + int2str(k[n]) + "]";
else
m_mult[rxn] += "*pow(c[" _ int2str(k[n]) + "],"+fp2str(order[n])+")";
if (stoich[n] == 1.0) {
m_is[k[n]] += " + r[" + int2str(rxn) + "]";
m_ds[k[n]] += " - r[" + int2str(rxn) + "]";
m_ir[rxn] += " + g[" + int2str(k[n]) + "]";
m_dr[rxn] += " - g[" + int2str(k[n]) + "]";
}
else {
s = fp2str(stoich[n]);
m_is[k[n]] += " + "+s+"*r[" + int2str(rxn) + "]";
m_ds[k[n]] += " - "+s+"*r[" + int2str(rxn) + "]";
m_ir[rxn] += " + "+s+"*g[" + int2str(k[n]) + "]";
m_dr[rxn] += " - "+s+"*g[" + int2str(k[n]) + "]";
}
}
}
string mult(int rxn) { return m_mult[rxn]; }
string incrSpec(int k) { return m_is[k]; }
string incrSpec(int k, string) { return m_is[k]; }
string decrSpec(int k) { return m_ds[k]; }
string incrRxn(int rxn) { return m_ir[rxn]; }
string decrRxn(int rxn) { return m_dr[rxn]; }
@ -632,3 +820,4 @@ namespace Cantera {
}
#endif

View file

@ -62,7 +62,7 @@ namespace Cantera {
}
void StoichSubstance::setParametersFromXML(const XML_Node& eosdata) {
eosdata.require("model","StoichSubstance");
eosdata._require("model","StoichSubstance");
doublereal rho = getFloat(eosdata, "density", "-");
setDensity(rho);
}

View file

@ -241,7 +241,7 @@ namespace Cantera {
void SurfPhase::
setParametersFromXML(const XML_Node& eosdata) {
eosdata.require("model","Surface");
eosdata._require("model","Surface");
doublereal n = getFloat(eosdata, "site_density", "-");
if (n <= 0.0)
throw CanteraError("SurfPhase::setParametersFromXML",
@ -271,7 +271,7 @@ namespace Cantera {
void EdgePhase::
setParametersFromXML(const XML_Node& eosdata) {
eosdata.require("model","Edge");
eosdata._require("model","Edge");
doublereal n = getFloat(eosdata, "site_density", "-");
if (n <= 0.0)
throw CanteraError("EdgePhase::setParametersFromXML",

View file

@ -129,6 +129,10 @@ namespace Cantera {
typedef ct::ctvector_fp vector_fp;
typedef ct::ctvector_int array_int;
typedef ct::ctvector_int vector_int;
//typedef std::valarray<double> array_fp;
//typedef std::valarray<double> vector_fp;
//typedef std::valarray<int> array_int;
//typedef std::valarray<int> vector_int;
typedef vector_int group_t;
typedef std::vector<group_t> grouplist_t;

View file

@ -12,7 +12,10 @@
#ifndef CT_CTLAPACK_H
#define CT_CTLAPACK_H
#ifdef DARWIN
#undef USE_CBLAS
#undef NO_FTN_STRING_LEN_AT_END
#endif
#include "ct_defs.h"
@ -30,6 +33,8 @@
#define _DGBTRF_ dgbtrf
#define _DGBTRS_ dgbtrs
#define _DSCAL_ dscal
#else
#define _DGEMV_ dgemv_
@ -41,20 +46,26 @@
#define _DGBTRF_ dgbtrf_
#define _DGBTRS_ dgbtrs_
#define _DSCAL_ dscal_
#endif
namespace ctlapack {
typedef enum {Transpose = 1, NoTranspose = 0} transpose_t;
typedef enum {ColMajor = 1, RowMajor = 0} storage_t;
}
const char no_yes[2] = {'N', 'T'};
//const CBLAS_ORDER cblasOrder[2] = { CblasRowMajor, CblasColMajor };
//const CBLAS_TRANSPOSE cblasTrans[2] = { CblasNoTrans, CblasTrans };
#ifdef USE_CBLAS
#include <Accelerate.h>
const CBLAS_ORDER cblasOrder[2] = { CblasRowMajor, CblasColMajor };
const CBLAS_TRANSPOSE cblasTrans[2] = { CblasNoTrans, CblasTrans };
#endif
//#ifdef DARWIN
//#include <Accelerate.h>
//#else
// C interfaces for Fortran Lapack routines
extern "C" {
@ -120,7 +131,11 @@ extern "C" {
doublereal *b, integer *ldb, integer *info);
#endif
int _DSCAL_(integer *n, doublereal *da, doublereal *dx, integer *incx);
void cblas_dscal(const int N, const double alpha, double *X, const int incX);
}
//#endif
namespace Cantera {
@ -130,10 +145,6 @@ namespace Cantera {
const doublereal* x, int incX, doublereal beta,
doublereal* y, int incY)
{
//#ifdef HAVE_INTEL_MKL
//cblas_dgemv(cblasOrder[storage], cblasTrans[trans], m, n, alpha,
// a, lda, x, incX, beta, y, incY);
//#else
#ifdef USE_CBLAS
cblas_dgemv(cblasOrder[storage], cblasTrans[trans], m, n, alpha,
a, lda, x, incX, beta, y, incY);
@ -141,14 +152,20 @@ namespace Cantera {
integer f_m = m, f_n = n, f_lda = lda, f_incX = incX, f_incY = incY;
doublereal f_alpha = alpha, f_beta = beta;
ftnlen trsize = 1;
#ifdef NO_FTN_STRING_LEN_AT_END
_DGEMV_(&no_yes[trans], &f_m, &f_n, &f_alpha, a,
&f_lda, x, &f_incX, &f_beta, y, &f_incY);
#else
#ifdef LAPACK_FTN_STRING_LEN_AT_END
_DGEMV_(&no_yes[trans], &f_m, &f_n, &f_alpha, a,
&f_lda, x, &f_incX, &f_beta, y, &f_incY, trsize);
#else
_DGEMV_(&no_yes[trans], trsize, &f_m, &f_n, &f_alpha, a,
&f_lda, x, &f_incX, &f_beta, y, &f_incY);
#endif
#endif
#endif
}
@ -176,6 +193,10 @@ namespace Cantera {
integer f_n = n, f_kl = kl, f_ku = ku, f_nrhs = nrhs, f_lda = lda,
f_ldb = ldb, f_info = info;
char tr = no_yes[trans];
#ifdef NO_FTN_STRING_LEN_AT_END
_DGBTRS_(&tr, &f_n, &f_kl, &f_ku, &f_nrhs, a, &f_lda, ipiv,
b, &f_ldb, &f_info);
#else
ftnlen trsize = 1;
#ifdef LAPACK_FTN_STRING_LEN_AT_END
_DGBTRS_(&tr, &f_n, &f_kl, &f_ku, &f_nrhs, a, &f_lda, ipiv,
@ -183,6 +204,7 @@ namespace Cantera {
#else
_DGBTRS_(&tr, trsize, &f_n, &f_kl, &f_ku, &f_nrhs, a, &f_lda, ipiv,
b, &f_ldb, &f_info);
#endif
#endif
info = f_info;
}
@ -205,6 +227,10 @@ namespace Cantera {
f_info = info;
char tr = no_yes[trans];
#ifdef NO_FTN_STRING_LEN_AT_END
_DGETRS_(&tr, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb,
&f_info);
#else
ftnlen trsize = 1;
#ifdef LAPACK_FTN_STRING_LEN_AT_END
_DGETRS_(&tr, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb,
@ -212,6 +238,7 @@ namespace Cantera {
#else
_DGETRS_(&tr, trsize, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb,
&f_info);
#endif
#endif
info = f_info;
}
@ -236,6 +263,13 @@ namespace Cantera {
rank = f_rank;
}
inline void ct_dscal(int n, doublereal da, doublereal* dx, int incx) {
//integer f_n = n, f_incx = incx;
//_DSCAL_(&f_n, &da, dx, &f_incx);
cblas_dscal(n, da, dx, incx);
}
}

View file

@ -728,8 +728,11 @@ namespace Cantera {
writelog("Logfile error."
"\n beginLogGroup: "+ __app->loggroups.back()+
"\n endLogGroup; "+title+"\n");
write_logfile("logerror");
__app->loggroups.clear();
__app->loglevels.clear();
}
if (__app->loggroups.size() == 1) {
else if (__app->loggroups.size() == 1) {
write_logfile(__app->loggroups.back()+"_log");
__app->loggroups.clear();
__app->loglevels.clear();

View file

@ -13,8 +13,8 @@
#define CT_DOMAIN1D_H
#include "../ctexceptions.h"
#include "../xml.h"
#include "../ctexceptions.h"
#include "refine.h"

View file

@ -15,11 +15,12 @@ INCDIR = ../../../build/include/cantera/kernel/oneD
INSTALL_TSC = ../../../bin/install_tsc
do_ranlib = @DO_RANLIB@
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
CXX_INCLUDES = -I..
CXX_INCLUDES = -I.. @CXX_INCLUDES@
# stirred reactors
OBJS = MultiJac.o MultiNewton.o newton_utils.o OneDim.o\
StFlow.o boundaries1D.o refine.o Sim1D.o
OBJS = oneD_files.o
#MultiJac.o MultiNewton.o newton_utils.o OneDim.o\
# StFlow.o boundaries1D.o refine.o Sim1D.o
ONED_H = Inlet1D.h MultiJac.h Sim1D.h StFlow.h \
Surf1D.h Domain1D.h MultiNewton.h OneDim.h \
Resid1D.h Solid1D.h refine.h

View file

@ -19,8 +19,8 @@
#pragma warning(disable:4503)
#endif
#include <vector>
#include <math.h>
//#include <vector>
//#include <math.h>
using namespace std;

View file

@ -13,7 +13,7 @@
namespace Cantera {
static void drawline() {
static void sim1D_drawline() {
string s(78,'.');
s += '\n';
writelog(s.c_str());
@ -250,7 +250,7 @@ namespace Cantera {
try {
if (loglevel > 0) {
drawline();
sim1D_drawline();
writelog("\nAttempt Newton solution of steady-state problem...");
}
newtonSolve(loglevel-1);
@ -276,7 +276,7 @@ namespace Cantera {
char buf[100];
if (loglevel > 0) {
writelog(" failure. \n\n");
drawline();
sim1D_drawline();
// }
//if (loglevel == 1)
writelog("Take "+int2str(nsteps)+

View file

@ -92,7 +92,7 @@ namespace Cantera {
}
static void drawline() {
static void st_drawline() {
writelog("\n-------------------------------------"
"------------------------------------------");
}
@ -890,14 +890,14 @@ namespace Cantera {
sprintf(buf, " Pressure: %10.4g Pa \n", m_press);
writelog(buf);
for (i = 0; i < nn; i++) {
drawline();
st_drawline();
sprintf(buf, "\n z ");
writelog(buf);
for (n = 0; n < 5; n++) {
sprintf(buf, " %10s ",componentName(i*5 + n).c_str());
writelog(buf);
}
drawline();
st_drawline();
for (j = 0; j < m_points; j++) {
sprintf(buf, "\n %10.4g ",m_z[j]);
writelog(buf);
@ -909,14 +909,14 @@ namespace Cantera {
writelog("\n");
}
int nrem = m_nv - 5*nn;
drawline();
st_drawline();
sprintf(buf, "\n z ");
writelog(buf);
for (n = 0; n < nrem; n++) {
sprintf(buf, " %10s ", componentName(nn*5 + n).c_str());
writelog(buf);
}
drawline();
st_drawline();
for (j = 0; j < m_points; j++) {
sprintf(buf, "\n %10.4g ",m_z[j]);
writelog(buf);

View file

@ -26,7 +26,7 @@ namespace Cantera {
* considers one domain.
*/
doublereal bound_step(const doublereal* x, const doublereal* step,
Domain1D& r, int loglevel=0) {
Domain1D& r, int loglevel) {
char buf[100];
int np = r.nPoints();

View file

@ -22,7 +22,7 @@ namespace Cantera {
return false;
}
static void drawline() {
static void r_drawline() {
string s(78,'#');
s += '\n';
writelog(s.c_str());
@ -204,7 +204,7 @@ namespace Cantera {
void Refiner::show() {
int nnew = static_cast<int>(m_loc.size());
if (nnew > 0) {
drawline();
r_drawline();
writelog(string("Refining grid in ") +
m_domain->id()+".\n"
+" New points inserted after grid points ");

View file

@ -182,6 +182,23 @@ namespace Cantera {
return logfile;
}
string wrapString(const string& s, int len) {
int nc = s.size();
int n, count=0;
string r;
for (n = 0; n < nc; n++) {
if (s[n] == '\n') count = 0;
else count++;
if (count > len && s[n] == ' ') {
r += "\n ";
count = 0;
}
r += s[n];
}
return r;
}
/*
* This routine strips off blanks and tabs (only leading and trailing
* characters) in 'str'. On return, it returns the number of
@ -303,5 +320,4 @@ namespace Cantera {
return rval;
}
}

View file

@ -36,6 +36,7 @@ namespace Cantera {
inline doublereal fpValue(string val) {
return atof(stripws(val).c_str());
}
string wrapString(const string& s, int len=70);
int stripLTWScstring(char str[]);
double atofCheck(const char *dptr);

View file

@ -18,13 +18,14 @@ do_ranlib = @DO_RANLIB@
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
# Transport Object Files
OBJS = TransportFactory.o MultiTransport.o MixTransport.o MMCollisionInt.o \
SolidTransport.o DustyGasTransport.o
OBJS = transport_files.o
#TransportFactory.o MultiTransport.o MixTransport.o MMCollisionInt.o \
# SolidTransport.o DustyGasTransport.o
TRAN_H = TransportFactory.h MultiTransport.h MixTransport.h \
MMCollisionInt.h SolidTransport.h DustyGasTransport.h \
TransportBase.h L_matrix.h FtnTransport.h TransportParams.h
CXX_INCLUDES = -I..
CXX_INCLUDES = -I.. @CXX_INCLUDES@
LIB = @buildlib@/libtransport.a
DEPENDS = $(OBJS:.o=.d)

View file

@ -82,9 +82,19 @@ namespace Cantera {
}
m_phi.resize(m_nsp, m_nsp, 0.0);
m_wratjk.resize(m_nsp, m_nsp, 0.0);
m_wratkj1.resize(m_nsp, m_nsp, 0.0);
int j, k;
for (j = 0; j < m_nsp; j++)
for (k = j; k < m_nsp; k++) {
m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
m_wratjk(k,j) = sqrt(m_wratjk(j,k));
m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
}
m_polytempvec.resize(5);
m_visc.resize(m_nsp);
m_sqvisc.resize(m_nsp);
m_cond.resize(m_nsp);
m_bdiff.resize(m_nsp, m_nsp);
@ -135,18 +145,15 @@ namespace Cantera {
if (m_viscmix_ok) return m_viscmix;
doublereal vismix = 0.0, denom;
int k, j;
doublereal vismix = 0.0;
int k;
// update m_visc and m_phi if necessary
if (!m_viscwt_ok) updateViscosity_T();
multiply(m_phi, m_molefracs.begin(), m_spwork.begin());
for (k = 0; k < m_nsp; k++) {
denom = 0.0;
for (j = 0; j < m_nsp; j++) {
denom += m_phi(k,j) * m_molefracs[j];
}
vismix += m_molefracs[k] * m_visc[k]/denom;
vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
}
m_viscmix = vismix;
return vismix;
@ -322,6 +329,7 @@ namespace Cantera {
m_logt = log(m_temp);
m_kbt = Boltzmann * m_temp;
m_sqrt_t = sqrt(m_temp);
m_t14 = sqrt(m_sqrt_t);
m_t32 = m_temp * m_sqrt_t;
m_sqrt_kbt = sqrt(Boltzmann*m_temp);
@ -438,13 +446,16 @@ namespace Cantera {
int k;
if (m_mode == CK_Mode) {
for (k = 0; k < m_nsp; k++) {
m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k]));
}
for (k = 0; k < m_nsp; k++) {
m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k]));
m_sqvisc[k] = sqrt(m_visc[k]);
}
}
else {
for (k = 0; k < m_nsp; k++) {
m_visc[k] = m_sqrt_t*dot5(m_polytempvec, m_visccoeffs[k]);
// the polynomial fit is done for sqrt(visc/sqrt(T))
m_sqvisc[k] = m_t14*dot5(m_polytempvec, m_visccoeffs[k]);
m_visc[k] = (m_sqvisc[k]*m_sqvisc[k]);
}
}
m_spvisc_ok = true;
@ -458,7 +469,7 @@ namespace Cantera {
* The flag m_visc_ok is set to true.
*/
void MixTransport::updateViscosity_T() {
doublereal vratiokj, wratiojk, rootwjk, factor1;
doublereal vratiokj, wratiojk, factor1;
if (!m_spvisc_ok) updateSpeciesViscosities();
@ -468,10 +479,12 @@ namespace Cantera {
for (k = j; k < m_nsp; k++) {
vratiokj = m_visc[k]/m_visc[j];
wratiojk = m_mw[j]/m_mw[k];
rootwjk = sqrt(wratiojk);
factor1 = 1.0 + sqrt(vratiokj * rootwjk);
// Note that m_wratjk(k,j) holds the square root of
// m_wratjk(j,k)!
factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
m_phi(k,j) = factor1*factor1 /
(SqrtEight * sqrt(1.0 + m_mw[k]/m_mw[j]));
(SqrtEight * m_wratkj1(j,k));
m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
}
}

View file

@ -114,6 +114,7 @@ namespace Cantera {
// property values
DenseMatrix m_bdiff;
vector_fp m_visc;
vector_fp m_sqvisc;
vector_fp m_cond;
array_fp m_molefracs;
@ -129,6 +130,7 @@ namespace Cantera {
DenseMatrix m_om22;
DenseMatrix m_phi; // viscosity weighting functions
DenseMatrix m_wratjk, m_wratkj1;
vector_fp m_zrot;
vector_fp m_crot;
@ -137,7 +139,7 @@ namespace Cantera {
vector_fp m_alpha;
vector_fp m_dipoleDiag;
doublereal m_temp, m_logt, m_kbt, m_t32;
doublereal m_temp, m_logt, m_kbt, m_t14, m_t32;
doublereal m_sqrt_kbt, m_sqrt_t;
vector_fp m_sqrt_eps_k;

View file

@ -25,7 +25,7 @@
#include "MultiTransport.h"
#include "../ctlapack.h"
#include "../../../ext/math/gmres.h"
//#include "../../../ext/math/gmres.h"
#include "../DenseMatrix.h"
#include "../polyfit.h"
@ -80,7 +80,7 @@ namespace Cantera {
/////////////////////////// constants //////////////////////////
const doublereal ThreeSixteenths = 3.0/16.0;
// const doublereal ThreeSixteenths = 3.0/16.0;
@ -182,10 +182,21 @@ namespace Cantera {
m_rotrelax.resize(m_nsp);
m_phi.resize(m_nsp, m_nsp, 0.0);
m_wratjk.resize(m_nsp, m_nsp, 0.0);
m_wratkj1.resize(m_nsp, m_nsp, 0.0);
int j, k;
for (j = 0; j < m_nsp; j++)
for (k = j; k < m_nsp; k++) {
m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
m_wratjk(k,j) = sqrt(m_wratjk(j,k));
m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
}
m_cinternal.resize(m_nsp);
m_polytempvec.resize(5);
m_visc.resize(m_nsp);
m_sqvisc.resize(m_nsp);
m_bdiff.resize(m_nsp, m_nsp);
//m_poly.resize(m_nsp);
@ -225,7 +236,7 @@ namespace Cantera {
// precompute and store log(epsilon_ij/k_B)
m_log_eps_k.resize(m_nsp, m_nsp);
int j;
// int j;
for (i = 0; i < m_nsp; i++) {
for (j = i; j < m_nsp; j++) {
m_log_eps_k(i,j) = log(tr.epsilon(i,j)/Boltzmann);
@ -239,7 +250,7 @@ namespace Cantera {
const doublereal sq298 = sqrt(298.0);
const doublereal kb298 = Boltzmann * 298.0;
m_sqrt_eps_k.resize(m_nsp);
int k;
//int k;
for (k = 0; k < m_nsp; k++) {
m_sqrt_eps_k[k] = sqrt(tr.eps[k]/Boltzmann);
m_frot_298[k] = Frot( tr.eps[k]/kb298,
@ -401,23 +412,26 @@ namespace Cantera {
// in m_a should provide a good starting guess, so convergence
// should be fast.
if (m_gmres) {
gmres(m_mgmres, 3*m_nsp, m_Lmatrix, m_b.begin(),
m_a.begin(), m_eps_gmres);
m_lmatrix_soln_ok = true;
m_l0000_ok = true; // L matrix not modified by GMRES
}
else {
//if (m_gmres) {
// gmres(m_mgmres, 3*m_nsp, m_Lmatrix, m_b.begin(),
// m_a.begin(), m_eps_gmres);
// m_lmatrix_soln_ok = true;
// m_l0000_ok = true; // L matrix not modified by GMRES
//}
//else {
copy(m_b.begin(), m_b.end(), m_a.begin());
int info = solve(m_Lmatrix, m_a.begin());
if (info != 0) {
try {
int info = solve(m_Lmatrix, m_a.begin());
}
catch (CanteraError) {
//if (info != 0) {
throw CanteraError("MultiTransport::solveLMatrixEquation",
"error in solving L matrix.");
}
m_lmatrix_soln_ok = true;
m_l0000_ok = false;
// L matrix is overwritten with LU decomposition
}
//}
m_lmatrix_soln_ok = true;
}
@ -757,6 +771,7 @@ namespace Cantera {
m_logt = log(m_temp);
m_kbt = Boltzmann * m_temp;
m_sqrt_t = sqrt(m_temp);
m_t14 = sqrt(m_sqrt_t);
m_t32 = m_temp * m_sqrt_t;
m_sqrt_kbt = sqrt(Boltzmann*m_temp);
@ -872,11 +887,15 @@ namespace Cantera {
if (m_mode == CK_Mode) {
for (k = 0; k < m_nsp; k++) {
m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k]));
m_sqvisc[k] = sqrt(m_visc[k]);
}
}
else {
for (k = 0; k < m_nsp; k++) {
m_visc[k] = m_sqrt_t*dot5(m_polytempvec, m_visccoeffs[k]);
//m_visc[k] = m_sqrt_t*dot5(m_polytempvec, m_visccoeffs[k]);
// the polynomial fit is done for sqrt(visc/sqrt(T))
m_sqvisc[k] = m_t14*dot5(m_polytempvec, m_visccoeffs[k]);
m_visc[k] = (m_sqvisc[k]*m_sqvisc[k]);
}
}
m_spvisc_ok = true;
@ -888,12 +907,11 @@ namespace Cantera {
void MultiTransport::updateViscosity_T() {
if (m_visc_tlast == m_thermo->temperature()) return;
_update_visc_T();
//m_thermo->update_T(m_update_visc_T);
m_visc_tlast = m_thermo->temperature();
}
void MultiTransport::_update_visc_T() {
doublereal vratiokj, wratiojk, rootwjk, factor1;
doublereal vratiokj, wratiojk, factor1;
updateSpeciesViscosities_T();
@ -903,10 +921,17 @@ namespace Cantera {
for (k = j; k < m_nsp; k++) {
vratiokj = m_visc[k]/m_visc[j];
wratiojk = m_mw[j]/m_mw[k];
rootwjk = sqrt(wratiojk);
factor1 = 1.0 + sqrt(vratiokj * rootwjk);
//rootwjk = sqrt(wratiojk);
//factor1 = 1.0 + sqrt(vratiokj * rootwjk);
//m_phi(k,j) = factor1*factor1 /
// (SqrtEight * sqrt(1.0 + m_mw[k]/m_mw[j]));
//m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
// Note that m_wratjk(k,j) holds the square root of
// m_wratjk(j,k)!
factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
m_phi(k,j) = factor1*factor1 /
(SqrtEight * sqrt(1.0 + m_mw[k]/m_mw[j]));
(SqrtEight * m_wratkj1(j,k));
m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
}
}

View file

@ -201,6 +201,7 @@ namespace Cantera {
// property values
DenseMatrix m_bdiff;
vector_fp m_visc;
vector_fp m_sqvisc;
array_fp m_molefracs;
@ -216,6 +217,7 @@ namespace Cantera {
DenseMatrix m_om22;
DenseMatrix m_phi; // viscosity weighting functions
DenseMatrix m_wratjk, m_wratkj1;
vector_fp m_zrot;
vector_fp m_crot;
@ -224,7 +226,7 @@ namespace Cantera {
vector_fp m_alpha;
vector_fp m_dipoleDiag;
doublereal m_temp, m_logt, m_kbt, m_t32;
doublereal m_temp, m_logt, m_kbt, m_t14, m_t32;
doublereal m_sqrt_kbt, m_sqrt_t;
vector_fp m_sqrt_eps_k;

View file

@ -808,7 +808,21 @@ namespace Cantera {
w2[n] = -1.0;
}
else {
spvisc[n] = visc/sqrt_T;
// the viscosity should be proportional
// approximately to sqrt(T); therefore,
// visc/sqrt(T) should have only a weak
// temperature dependence. And since the mixture
// rule requires the square root of the
// pure-species viscosity, fit the square root of
// (visc/sqrt(T)) to avoid having to compute
// square roots in the mixture rule.
spvisc[n] = sqrt(visc/sqrt_T);
// the pure-species conductivity scales
// approximately with sqrt(T). Unlike the
// viscosity, there is no reason here to fit the
// square root, since a different mixture rule is
// used.
spcond[n] = cond/sqrt_T;
w[n] = 1.0/(spvisc[n]*spvisc[n]);
w2[n] = 1.0/(spcond[n]*spcond[n]);
@ -827,8 +841,8 @@ namespace Cantera {
}
else {
sqrt_T = exp(0.5*tlog[n]);
val = sqrt_T * spvisc[n];
fit = sqrt_T * poly4(tlog[n], c.begin());
val = sqrt_T * pow(spvisc[n],2);
fit = sqrt_T * pow(poly4(tlog[n], c.begin()),2);
}
err = fit - val;
relerr = err/val;

View file

@ -11,6 +11,10 @@
#include "ct_defs.h"
#ifdef DARWINNN
#include <Accelerate.h>
#endif
namespace Cantera {
/**
@ -248,6 +252,14 @@ namespace Cantera {
return sum;
}
inline void scale(int N, double alpha, double* x) {
//#ifdef DARWINNNN
//cblas_dscal(N, alpha, x, 1);
//#else
scale(x, x+N, x, alpha);
//#endif
}
}

View file

@ -721,7 +721,7 @@ namespace Cantera {
}
void XML_Node::require(string a, string v) const {
void XML_Node::_require(string a, string v) const {
if (hasAttrib(a)) {
if (attrib(a) == v) return;
}

View file

@ -13,8 +13,8 @@
// Copyright 2001 California Institute of Technology
#ifndef CT_XML
#define CT_XML
#ifndef CT_XML_H
#define CT_XML_H
#include <string>
#include <vector>
@ -131,7 +131,9 @@ namespace Cantera {
int nChildren() const { return m_nchildren; }
void build(istream& f);
void require(string a, string v) const;
void _require(string a, string v) const;
/**
* This routine carries out a search for an XML node based
* on both the xml element name and the attribute ID.

16
config/configure vendored
View file

@ -271,7 +271,7 @@ PACKAGE_STRING=
PACKAGE_BUGREPORT=
ac_unique_file="Cantera.README"
ac_subst_vars='SHELL PATH_SEPARATOR PACKAGE_NAME PACKAGE_TARNAME PACKAGE_VERSION PACKAGE_STRING PACKAGE_BUGREPORT exec_prefix prefix program_transform_name bindir sbindir libexecdir datadir sysconfdir sharedstatedir localstatedir libdir includedir oldincludedir infodir mandir build_alias host_alias target_alias DEFS ECHO_C ECHO_N ECHO_T LIBS CVF_LIBDIR USE_CLIB_DLL local_inst local_python_inst python_prefix python_win_prefix ctversion homedir ct_libdir ct_bindir ct_incdir ct_incroot ct_datadir ct_demodir ct_templdir ct_tutdir ct_docdir ct_dir ct_mandir COMPACT_INSTALL build build_cpu build_vendor build_os host host_cpu host_vendor host_os target target_cpu target_vendor target_os username ctroot buildinc buildlib buildbin MAKE ARCHIVE DO_RANLIB RANLIB SOEXT SHARED PIC LCXX_FLAGS LCXX_END_LIBS USERDIR INCL_USER_CODE phase_object_files phase_header_files KERNEL KERNEL_OBJ BUILD_CK LIB_DIR build_lapack build_blas BLAS_LAPACK_LIBS BLAS_LAPACK_DIR build_with_f2c LOCAL_LIB_DIRS LOCAL_LIBS CT_SHARED_LIB F77FLAGS PYTHON_CMD BUILD_PYTHON NUMARRAY_INC_DIR NUMARRAY_HOME CANTERA_PYTHON_HOME CVSTAG MATLAB_CMD BUILD_MATLAB BUILD_CLIB export_name INSTALL_PROGRAM INSTALL_SCRIPT INSTALL_DATA CXX CXXFLAGS LDFLAGS CPPFLAGS ac_ct_CXX EXEEXT OBJEXT CC CFLAGS ac_ct_CC F77 FFLAGS ac_ct_F77 FLIBS F90 BUILD_F90 F90FLAGS F90BUILDFLAGS precompile_headers CXX_DEPENDS OS_IS_DARWIN OS_IS_WIN OS_IS_CYGWIN SHARED_CTLIB mex_ext F77_EXT CXX_EXT OBJ_EXT EXE_EXT local_math_libs math_libs SO LDSHARED EXTRA_LINK LIBOBJS LTLIBOBJS'
ac_subst_vars='SHELL PATH_SEPARATOR PACKAGE_NAME PACKAGE_TARNAME PACKAGE_VERSION PACKAGE_STRING PACKAGE_BUGREPORT exec_prefix prefix program_transform_name bindir sbindir libexecdir datadir sysconfdir sharedstatedir localstatedir libdir includedir oldincludedir infodir mandir build_alias host_alias target_alias DEFS ECHO_C ECHO_N ECHO_T LIBS CVF_LIBDIR USE_CLIB_DLL local_inst local_python_inst python_prefix python_win_prefix ctversion homedir ct_libdir ct_bindir ct_incdir ct_incroot ct_datadir ct_demodir ct_templdir ct_tutdir ct_docdir ct_dir ct_mandir COMPACT_INSTALL build build_cpu build_vendor build_os host host_cpu host_vendor host_os target target_cpu target_vendor target_os username ctroot buildinc buildlib buildbin MAKE ARCHIVE DO_RANLIB RANLIB SOEXT SHARED PIC LCXX_FLAGS LCXX_END_LIBS CXX_INCLUDES USERDIR INCL_USER_CODE phase_object_files phase_header_files KERNEL KERNEL_OBJ BUILD_CK LIB_DIR build_lapack build_blas BLAS_LAPACK_LIBS BLAS_LAPACK_DIR build_with_f2c LOCAL_LIB_DIRS LOCAL_LIBS CT_SHARED_LIB F77FLAGS PYTHON_CMD BUILD_PYTHON NUMARRAY_INC_DIR NUMARRAY_HOME CANTERA_PYTHON_HOME CVSTAG MATLAB_CMD BUILD_MATLAB BUILD_CLIB export_name INSTALL_PROGRAM INSTALL_SCRIPT INSTALL_DATA CXX CXXFLAGS LDFLAGS CPPFLAGS ac_ct_CXX EXEEXT OBJEXT CC CFLAGS ac_ct_CC F77 FFLAGS ac_ct_F77 FLIBS F90 BUILD_F90 F90FLAGS F90BUILDFLAGS precompile_headers CXX_DEPENDS OS_IS_DARWIN OS_IS_WIN OS_IS_CYGWIN SHARED_CTLIB mex_ext F77_EXT CXX_EXT OBJ_EXT EXE_EXT local_math_libs math_libs SO LDSHARED EXTRA_LINK LIBOBJS LTLIBOBJS'
ac_subst_files=''
# Initialize some variables set by options.
@ -1272,7 +1272,10 @@ EXTRA_LINK=""
mex_ext=mexglx
case $ac_sys_system in
Darwin*) OS_IS_DARWIN=1; EXTRA_LINK="-framework Accelerate"; mex_ext=mexmac;;
Darwin*) OS_IS_DARWIN=1;
EXTRA_LINK="-framework Accelerate";
CXX_INCLUDES="$CXX_INCLUDES -I/System/Library/Frameworks/Accelerate.framework/Headers"
mex_ext=mexmac;;
CYGWIN*) OS_IS_CYGWIN=1; mex_ext=dll;;
esac
@ -1571,6 +1574,8 @@ if test -z "$PIC"; then PIC='-fPIC'; fi
#if test -z "$LCXX_END_LIBS"; then LCXX_END_LIBS='-lm'; fi
#########################################################
# User Code
#########################################################
@ -3433,7 +3438,7 @@ fi
# Provide some information about the compiler.
echo "$as_me:3436:" \
echo "$as_me:3441:" \
"checking for Fortran 77 compiler version" >&5
ac_compiler=`set X $ac_compile; echo $2`
{ (eval echo "$as_me:$LINENO: \"$ac_compiler --version </dev/null >&5\"") >&5
@ -3610,7 +3615,7 @@ _ACEOF
# flags.
ac_save_FFLAGS=$FFLAGS
FFLAGS="$FFLAGS $ac_verb"
(eval echo $as_me:3613: \"$ac_link\") >&5
(eval echo $as_me:3618: \"$ac_link\") >&5
ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'`
echo "$ac_f77_v_output" >&5
FFLAGS=$ac_save_FFLAGS
@ -3690,7 +3695,7 @@ _ACEOF
# flags.
ac_save_FFLAGS=$FFLAGS
FFLAGS="$FFLAGS $ac_cv_prog_f77_v"
(eval echo $as_me:3693: \"$ac_link\") >&5
(eval echo $as_me:3698: \"$ac_link\") >&5
ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'`
echo "$ac_f77_v_output" >&5
FFLAGS=$ac_save_FFLAGS
@ -4762,6 +4767,7 @@ s,@SHARED@,$SHARED,;t t
s,@PIC@,$PIC,;t t
s,@LCXX_FLAGS@,$LCXX_FLAGS,;t t
s,@LCXX_END_LIBS@,$LCXX_END_LIBS,;t t
s,@CXX_INCLUDES@,$CXX_INCLUDES,;t t
s,@USERDIR@,$USERDIR,;t t
s,@INCL_USER_CODE@,$INCL_USER_CODE,;t t
s,@phase_object_files@,$phase_object_files,;t t

View file

@ -28,7 +28,10 @@ EXTRA_LINK=""
mex_ext=mexglx
case $ac_sys_system in
Darwin*) OS_IS_DARWIN=1; EXTRA_LINK="-framework Accelerate"; mex_ext=mexmac;;
Darwin*) OS_IS_DARWIN=1;
EXTRA_LINK="-framework Accelerate";
CXX_INCLUDES="$CXX_INCLUDES -I/System/Library/Frameworks/Accelerate.framework/Headers"
mex_ext=mexmac;;
CYGWIN*) OS_IS_CYGWIN=1; mex_ext=dll;;
esac
@ -228,6 +231,8 @@ AC_SUBST(LCXX_FLAGS)
#if test -z "$LCXX_END_LIBS"; then LCXX_END_LIBS='-lm'; fi
AC_SUBST(LCXX_END_LIBS)
AC_SUBST(CXX_INCLUDES)
#########################################################
# User Code
#########################################################

6
configure vendored
View file

@ -139,7 +139,7 @@ F90=${F90:="default"}
# these compilers will be added automatically, and you do not need to
# specify them here. Otherwise, add any required compiler-specific
# flags here.
F90FLAGS=${F90FLAGS:='-O3'}
F90FLAGS=${F90FLAGS:='-O3 -g'}
#----------------------------------------------------------------------
@ -261,7 +261,7 @@ CXX=${CXX:=g++}
CC=${CC:=gcc}
# C++ compiler flags
CXXFLAGS=${CXXFLAGS:="-O0 -g -Wall"}
CXXFLAGS=${CXXFLAGS:="-O3 -g -Wall"}
# the C++ flags required for linking. Uncomment if additional flags
# need to be passed to the linker.
@ -305,7 +305,7 @@ F77=${F77:=g77}
# Fortran 77 compiler flags. Note that the Fortran compiler flags must be set
# to produce object code compatible with the C/C++ compiler you are using.
FFLAGS=${FFLAGS:='-O0 -g'}
FFLAGS=${FFLAGS:='-O3 -g'}
# the additional Fortran flags required for linking, if any. Leave commented
# out if no additional flags are required.

View file

@ -37,7 +37,7 @@ FORT_LIBS = @FLIBS@
CXX = @CXX@
# C++ compile flags
CXX_FLAGS = @CXXFLAGS@
CXX_FLAGS = @CXXFLAGS@ @CXX_INCLUDES@
# external libraries
EXT_LIBS = @LOCAL_LIBS@ -lctcxx

View file

@ -76,8 +76,7 @@ extern "C" {
******************************************************************/
typedef double real;
/* typedef long integer; dgg */
//typedef long int integer;
typedef int integer;
#define LLNL_FLOAT 0

View file

@ -203,8 +203,8 @@ void dswap_( const int *n, double *x, const int *incx, double *y,
const int *incy );
// x <= a*x
extern "C"
void dscal_( const int *n, const double *alpha, double *x, const int *incx );
//extern "C"
//void dscal_( const int *n, const double *alpha, double *x, const int *incx );
// y <= x

View file

@ -9,7 +9,7 @@ LCXX_FLAGS = -L$(LIBDIR) @CXXFLAGS@
LOCAL_LIBS = -lcantera -ltpx -lctcxx
#LOCAL_LIBS = -lcantera @math_libs@ @BLAS_LAPACK_LIBS@ -lctcxx
LCXX_END_LIBS = @LCXX_END_LIBS@
LCXX_END_LIBS = @LCXX_END_LIBS@ @EXTRA_LINK@
OBJS = ck2cti.o cti2ctml.o fixtext.o

View file

@ -30,7 +30,7 @@ LINK_OPTIONS = @LCXX_FLAGS@ @EXTRA_LINK@
CXX = @CXX@
# C++ compile flags
CXX_FLAGS = @CXXFLAGS@
CXX_FLAGS = @CXXFLAGS@ @CXX_INCLUDES@
# external libraries
EXT_LIBS = @LOCAL_LIBS@ -lctcxx