diff --git a/Cantera/clib/src/ctmultiphase.cpp b/Cantera/clib/src/ctmultiphase.cpp index e858b15d1..929ee365e 100644 --- a/Cantera/clib/src/ctmultiphase.cpp +++ b/Cantera/clib/src/ctmultiphase.cpp @@ -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); + } + } diff --git a/Cantera/clib/src/ctmultiphase.h b/Cantera/clib/src/ctmultiphase.h index f85ea8fe1..614642b3a 100644 --- a/Cantera/clib/src/ctmultiphase.h +++ b/Cantera/clib/src/ctmultiphase.h @@ -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 diff --git a/Cantera/clib/src/ctonedim.cpp b/Cantera/clib/src/ctonedim.cpp index f633d9986..0adda95fe 100644 --- a/Cantera/clib/src/ctonedim.cpp +++ b/Cantera/clib/src/ctonedim.cpp @@ -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; diff --git a/Cantera/clib/src/ctonedim.h b/Cantera/clib/src/ctonedim.h index 6a6ebb3e8..c5e1521b5 100644 --- a/Cantera/clib/src/ctonedim.h +++ b/Cantera/clib/src/ctonedim.h @@ -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); diff --git a/Cantera/cxx/demos/Makefile.in b/Cantera/cxx/demos/Makefile.in index 1d4aab328..e71afe5cf 100644 --- a/Cantera/cxx/demos/Makefile.in +++ b/Cantera/cxx/demos/Makefile.in @@ -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 diff --git a/Cantera/python/Cantera/mixture.py b/Cantera/python/Cantera/mixture.py index 66556e39b..46931592b 100644 --- a/Cantera/python/Cantera/mixture.py +++ b/Cantera/python/Cantera/mixture.py @@ -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 diff --git a/Cantera/python/src/ctmultiphase_methods.cpp b/Cantera/python/src/ctmultiphase_methods.cpp index 1330d88b1..9c6c933c8 100644 --- a/Cantera/python/src/ctmultiphase_methods.cpp +++ b/Cantera/python/src/ctmultiphase_methods.cpp @@ -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) { diff --git a/Cantera/python/src/ctonedim_methods.cpp b/Cantera/python/src/ctonedim_methods.cpp index daf20594e..42ffb0b49 100644 --- a/Cantera/python/src/ctonedim_methods.cpp +++ b/Cantera/python/src/ctonedim_methods.cpp @@ -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); diff --git a/Cantera/python/src/methods.h b/Cantera/python/src/methods.h index 40394a2fd..efeececea 100644 --- a/Cantera/python/src/methods.h +++ b/Cantera/python/src/methods.h @@ -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}, diff --git a/Cantera/src/ConstDensityThermo.cpp b/Cantera/src/ConstDensityThermo.cpp index 74416eb67..57f33a2c8 100755 --- a/Cantera/src/ConstDensityThermo.cpp +++ b/Cantera/src/ConstDensityThermo.cpp @@ -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); } diff --git a/Cantera/src/DenseMatrix.cpp b/Cantera/src/DenseMatrix.cpp index 6b65ebe44..3a1a1f477 100755 --- a/Cantera/src/DenseMatrix.cpp +++ b/Cantera/src/DenseMatrix.cpp @@ -108,7 +108,7 @@ namespace Cantera { if (info != 0) throw CanteraError("DenseMatrix::leaseSquares", "DGELSS returned INFO = "+int2str(info)); - return 0; + return 0; } #endif diff --git a/Cantera/src/GasKinetics.cpp b/Cantera/src/GasKinetics.cpp index 7c813a869..5303f242d 100755 --- a/Cantera/src/GasKinetics.cpp +++ b/Cantera/src/GasKinetics.cpp @@ -21,6 +21,8 @@ #include "ThirdBodyMgr.h" #include "RateCoeffMgr.h" +//#include "../user/grirxnstoich.h" + #include 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; } } diff --git a/Cantera/src/GasKinetics.h b/Cantera/src/GasKinetics.h index d48789e22..e1ecc7908 100755 --- a/Cantera/src/GasKinetics.h +++ b/Cantera/src/GasKinetics.h @@ -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 m_fwdOrder; diff --git a/Cantera/src/LatticeSolidPhase.cpp b/Cantera/src/LatticeSolidPhase.cpp index 0a836f903..46768fea9 100644 --- a/Cantera/src/LatticeSolidPhase.cpp +++ b/Cantera/src/LatticeSolidPhase.cpp @@ -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 lattices; la.getChildren("Lattice",lattices); diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in index 98e039d1a..24f9d18df 100755 --- a/Cantera/src/Makefile.in +++ b/Cantera/src/Makefile.in @@ -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) diff --git a/Cantera/src/MetalPhase.h b/Cantera/src/MetalPhase.h index 9186b2a3f..9e7dd3a43 100644 --- a/Cantera/src/MetalPhase.h +++ b/Cantera/src/MetalPhase.h @@ -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); } diff --git a/Cantera/src/MultiPhase.cpp b/Cantera/src/MultiPhase.cpp index 373e03c30..45ba0f96e 100644 --- a/Cantera/src/MultiPhase.cpp +++ b/Cantera/src/MultiPhase.cpp @@ -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", diff --git a/Cantera/src/MultiPhase.h b/Cantera/src/MultiPhase.h index fae29569d..44e4148bf 100644 --- a/Cantera/src/MultiPhase.h +++ b/Cantera/src/MultiPhase.h @@ -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(); diff --git a/Cantera/src/MultiPhaseEquil.cpp b/Cantera/src/MultiPhaseEquil.cpp index e99dbd581..36627ad23 100644 --- a/Cantera/src/MultiPhaseEquil.cpp +++ b/Cantera/src/MultiPhaseEquil.cpp @@ -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; } diff --git a/Cantera/src/PureFluidPhase.cpp b/Cantera/src/PureFluidPhase.cpp index 5c03c3680..e9d0a7c7d 100644 --- a/Cantera/src/PureFluidPhase.cpp +++ b/Cantera/src/PureFluidPhase.cpp @@ -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", diff --git a/Cantera/src/ReactionStoichMgr.cpp b/Cantera/src/ReactionStoichMgr.cpp index 68755946c..0fb2506af 100644 --- a/Cantera/src/ReactionStoichMgr.cpp +++ b/Cantera/src/ReactionStoichMgr.cpp @@ -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 out; + m_revproducts->writeIncrementSpecies("rf",out); + m_irrevproducts->writeIncrementSpecies("rf",out); + m_reactants->writeIncrementSpecies("rb",out); + map::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 out; + m_revproducts->writeIncrementSpecies("rb",out); + m_reactants->writeIncrementSpecies("rf",out); + map::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 out; + m_revproducts->writeIncrementSpecies("r",out); + m_irrevproducts->writeIncrementSpecies("r",out); + m_reactants->writeDecrementSpecies("r",out); + map::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 out; + m_reactants->writeMultiply("c",out); + map::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 out; + m_revproducts->writeMultiply("c",out); + map::iterator b; + for (b = out.begin(); b != out.end(); ++b) { + string rhs = b->second; + f << " r[" << b->first << "] *= " << rhs << ";" << endl; + } + f << " }" << endl << endl << endl; + } } diff --git a/Cantera/src/ReactionStoichMgr.h b/Cantera/src/ReactionStoichMgr.h index 3ad208754..e30600ed0 100644 --- a/Cantera/src/ReactionStoichMgr.h +++ b/Cantera/src/ReactionStoichMgr.h @@ -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 }; } diff --git a/Cantera/src/State.cpp b/Cantera/src/State.cpp index 574602c9e..721985d3b 100644 --- a/Cantera/src/State.cpp +++ b/Cantera/src/State.cpp @@ -19,7 +19,9 @@ #include "ctexceptions.h" #include "stringUtils.h" #include "State.h" - +#ifdef DARWIN +#include +#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 { } } - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/Cantera/src/StoichManager.h b/Cantera/src/StoichManager.h index a71ee36cd..e4054325f 100755 --- a/Cantera/src/StoichManager.h +++ b/Cantera/src/StoichManager.h @@ -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& out) { + out[m_rxn] = fmt(r, m_ic0); + } + + void writeIncrementReaction(string r, map& out) { + out[m_rxn] += " + "+fmt(r, m_ic0); + } + void writeDecrementReaction(string r, map& out) { + out[m_rxn] += " - "+fmt(r, m_ic0); + } + + void writeIncrementSpecies(string r, map& out) { + out[m_ic0] += " + "+fmt(r, m_rxn); + } + void writeDecrementSpecies(string r, map& 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& out) { + out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1); + } + void writeIncrementReaction(string r, map& out) { + out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1); + } + void writeDecrementReaction(string r, map& 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& out) { + string s = " + "+fmt(r, m_rxn); + out[m_ic0] += s; + out[m_ic1] += s; + } + void writeDecrementSpecies(string r, map& 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& out) { + out[m_rxn] = fmt(r, m_ic0) + " * " + fmt(r, m_ic1) + " * " + fmt(r, m_ic2); + } + void writeIncrementReaction(string r, map& out) { + out[m_rxn] += " + "+fmt(r, m_ic0)+" + "+fmt(r, m_ic1)+" + "+fmt(r, m_ic2); + } + void writeDecrementReaction(string r, map& out) { + out[m_rxn] += " - "+fmt(r, m_ic0)+" - "+fmt(r, m_ic1)+" - "+fmt(r, m_ic2); + } + void writeIncrementSpecies(string r, map& out) { + string s = " + "+fmt(r, m_rxn); + out[m_ic0] += s; + out[m_ic1] += s; + out[m_ic2] += s; + } + void writeDecrementSpecies(string r, map& 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& 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& 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& 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& 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& 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 + inline static void _writeIncrementSpecies(_InputIter __begin, _InputIter __end, string r, + map& out) { + for (; __begin != __end; ++__begin) __begin->writeIncrementSpecies(r, out); + } + + template + inline static void _writeDecrementSpecies(_InputIter __begin, _InputIter __end, string r, + map& out) { + for (; __begin != __end; ++__begin) __begin->writeDecrementSpecies(r, out); + } + + template + inline static void _writeIncrementReaction(_InputIter __begin, _InputIter __end, string r, + map& out) { + for (; __begin != __end; ++__begin) __begin->writeIncrementReaction(r, out); + } + + template + inline static void _writeDecrementReaction(_InputIter __begin, _InputIter __end, string r, + map& out) { + for (; __begin != __end; ++__begin) __begin->writeDecrementReaction(r, out); + } + + template + inline static void _writeMultiply(_InputIter __begin, _InputIter __end, string r, + map& 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& 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& 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& 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& 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& 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 m_c1_list; @@ -590,7 +760,7 @@ namespace Cantera { map 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 + diff --git a/Cantera/src/StoichSubstance.cpp b/Cantera/src/StoichSubstance.cpp index 620d1d977..c2eff8024 100644 --- a/Cantera/src/StoichSubstance.cpp +++ b/Cantera/src/StoichSubstance.cpp @@ -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); } diff --git a/Cantera/src/SurfPhase.cpp b/Cantera/src/SurfPhase.cpp index 3f94bf40b..8242b8e3d 100644 --- a/Cantera/src/SurfPhase.cpp +++ b/Cantera/src/SurfPhase.cpp @@ -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", diff --git a/Cantera/src/ct_defs.h b/Cantera/src/ct_defs.h index 4379fdcb4..16c37c00c 100755 --- a/Cantera/src/ct_defs.h +++ b/Cantera/src/ct_defs.h @@ -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 array_fp; + //typedef std::valarray vector_fp; + //typedef std::valarray array_int; + //typedef std::valarray vector_int; typedef vector_int group_t; typedef std::vector grouplist_t; diff --git a/Cantera/src/ctlapack.h b/Cantera/src/ctlapack.h index 792cdfc6e..436fdb9c0 100755 --- a/Cantera/src/ctlapack.h +++ b/Cantera/src/ctlapack.h @@ -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 +const CBLAS_ORDER cblasOrder[2] = { CblasRowMajor, CblasColMajor }; +const CBLAS_TRANSPOSE cblasTrans[2] = { CblasNoTrans, CblasTrans }; +#endif +//#ifdef DARWIN +//#include +//#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); + + } + } diff --git a/Cantera/src/misc.cpp b/Cantera/src/misc.cpp index 86663f6ce..e84a6955c 100755 --- a/Cantera/src/misc.cpp +++ b/Cantera/src/misc.cpp @@ -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(); diff --git a/Cantera/src/oneD/Domain1D.h b/Cantera/src/oneD/Domain1D.h index 33122fdf4..c7388f68f 100644 --- a/Cantera/src/oneD/Domain1D.h +++ b/Cantera/src/oneD/Domain1D.h @@ -13,8 +13,8 @@ #define CT_DOMAIN1D_H -#include "../ctexceptions.h" #include "../xml.h" +#include "../ctexceptions.h" #include "refine.h" diff --git a/Cantera/src/oneD/Makefile.in b/Cantera/src/oneD/Makefile.in index 93a9f4ebd..01f3d6b91 100644 --- a/Cantera/src/oneD/Makefile.in +++ b/Cantera/src/oneD/Makefile.in @@ -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 diff --git a/Cantera/src/oneD/MultiJac.cpp b/Cantera/src/oneD/MultiJac.cpp index f65d1f8f1..3d69b6be0 100644 --- a/Cantera/src/oneD/MultiJac.cpp +++ b/Cantera/src/oneD/MultiJac.cpp @@ -19,8 +19,8 @@ #pragma warning(disable:4503) #endif -#include -#include +//#include +//#include using namespace std; diff --git a/Cantera/src/oneD/Sim1D.cpp b/Cantera/src/oneD/Sim1D.cpp index b1f8619c1..b91662e28 100644 --- a/Cantera/src/oneD/Sim1D.cpp +++ b/Cantera/src/oneD/Sim1D.cpp @@ -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)+ diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp index 60d9766b9..8593a94cf 100644 --- a/Cantera/src/oneD/StFlow.cpp +++ b/Cantera/src/oneD/StFlow.cpp @@ -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); diff --git a/Cantera/src/oneD/newton_utils.cpp b/Cantera/src/oneD/newton_utils.cpp index 4690d9177..93042c75d 100644 --- a/Cantera/src/oneD/newton_utils.cpp +++ b/Cantera/src/oneD/newton_utils.cpp @@ -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(); diff --git a/Cantera/src/oneD/refine.cpp b/Cantera/src/oneD/refine.cpp index 23c610747..5626e93a3 100644 --- a/Cantera/src/oneD/refine.cpp +++ b/Cantera/src/oneD/refine.cpp @@ -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(m_loc.size()); if (nnew > 0) { - drawline(); + r_drawline(); writelog(string("Refining grid in ") + m_domain->id()+".\n" +" New points inserted after grid points "); diff --git a/Cantera/src/stringUtils.cpp b/Cantera/src/stringUtils.cpp index b5caaa10c..6077d6f0f 100755 --- a/Cantera/src/stringUtils.cpp +++ b/Cantera/src/stringUtils.cpp @@ -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; } - } diff --git a/Cantera/src/stringUtils.h b/Cantera/src/stringUtils.h index 4bd57cbc8..ebe4a0c7a 100755 --- a/Cantera/src/stringUtils.h +++ b/Cantera/src/stringUtils.h @@ -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); diff --git a/Cantera/src/transport/Makefile.in b/Cantera/src/transport/Makefile.in index 4480868dc..e5b407348 100644 --- a/Cantera/src/transport/Makefile.in +++ b/Cantera/src/transport/Makefile.in @@ -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) diff --git a/Cantera/src/transport/MixTransport.cpp b/Cantera/src/transport/MixTransport.cpp index da481caaa..bfd4ffb14 100755 --- a/Cantera/src/transport/MixTransport.cpp +++ b/Cantera/src/transport/MixTransport.cpp @@ -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); } } diff --git a/Cantera/src/transport/MixTransport.h b/Cantera/src/transport/MixTransport.h index f0237b2e2..d98b844cb 100755 --- a/Cantera/src/transport/MixTransport.h +++ b/Cantera/src/transport/MixTransport.h @@ -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; diff --git a/Cantera/src/transport/MultiTransport.cpp b/Cantera/src/transport/MultiTransport.cpp index 5b3988e7d..2702e834b 100755 --- a/Cantera/src/transport/MultiTransport.cpp +++ b/Cantera/src/transport/MultiTransport.cpp @@ -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); } } diff --git a/Cantera/src/transport/MultiTransport.h b/Cantera/src/transport/MultiTransport.h index d176a68b5..212ed66d0 100755 --- a/Cantera/src/transport/MultiTransport.h +++ b/Cantera/src/transport/MultiTransport.h @@ -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; diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index e8f7996e2..8bb56c138 100755 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -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; diff --git a/Cantera/src/utilities.h b/Cantera/src/utilities.h index 9676430d3..1eea7258a 100755 --- a/Cantera/src/utilities.h +++ b/Cantera/src/utilities.h @@ -11,6 +11,10 @@ #include "ct_defs.h" +#ifdef DARWINNN +#include +#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 + } + } diff --git a/Cantera/src/xml.cpp b/Cantera/src/xml.cpp index 70207788c..a6bbecad7 100755 --- a/Cantera/src/xml.cpp +++ b/Cantera/src/xml.cpp @@ -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; } diff --git a/Cantera/src/xml.h b/Cantera/src/xml.h index 0e7e57e13..54dc1cdef 100755 --- a/Cantera/src/xml.h +++ b/Cantera/src/xml.h @@ -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 #include @@ -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. diff --git a/config/configure b/config/configure index f1c2e18bf..865f8ab5d 100755 --- a/config/configure +++ b/config/configure @@ -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 &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 diff --git a/config/configure.in b/config/configure.in index a07a6a334..6d3f373cd 100755 --- a/config/configure.in +++ b/config/configure.in @@ -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 ######################################################### diff --git a/configure b/configure index e96b57e18..609dce4ed 100755 --- a/configure +++ b/configure @@ -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. diff --git a/examples/cxx/Makefile.in b/examples/cxx/Makefile.in index eab28b222..bf04ed4a9 100755 --- a/examples/cxx/Makefile.in +++ b/examples/cxx/Makefile.in @@ -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 diff --git a/ext/cvode/include/llnltyps.h b/ext/cvode/include/llnltyps.h index 5695625ed..20764e9cc 100755 --- a/ext/cvode/include/llnltyps.h +++ b/ext/cvode/include/llnltyps.h @@ -76,8 +76,7 @@ extern "C" { ******************************************************************/ typedef double real; - -/* typedef long integer; dgg */ + //typedef long int integer; typedef int integer; #define LLNL_FLOAT 0 diff --git a/ext/math/cblas.h b/ext/math/cblas.h index 8ca35aa1a..204698e97 100755 --- a/ext/math/cblas.h +++ b/ext/math/cblas.h @@ -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 diff --git a/tools/src/Makefile.in b/tools/src/Makefile.in index 3eadd1913..26d2a99d1 100755 --- a/tools/src/Makefile.in +++ b/tools/src/Makefile.in @@ -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 diff --git a/tools/templates/cxx/demo.mak.in b/tools/templates/cxx/demo.mak.in index be2882225..7b5ff4da7 100644 --- a/tools/templates/cxx/demo.mak.in +++ b/tools/templates/cxx/demo.mak.in @@ -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