diff --git a/Cantera/src/IdealGasThermo.cpp b/Cantera/src/IdealGasThermo.cpp deleted file mode 100755 index 942f88135..000000000 --- a/Cantera/src/IdealGasThermo.cpp +++ /dev/null @@ -1,108 +0,0 @@ -/** - * - * @file IdealGasThermo.cpp - * - */ - -#ifdef WIN32 -#pragma warning(disable:4786) -#pragma warning(disable:4503) -#endif - -#include "ct_defs.h" -#include "mix_defs.h" -#include "IdealGasThermo.h" -#include "SpeciesThermo.h" - -namespace Cantera { - - void IdealGasThermo::getChemPotentials(doublereal* mu) const { - doublereal logp = log(pressure()/m_spthermo->refPressure()); - doublereal xx; - doublereal rt = m_s->temperature() * GasConstant; - const array_fp& g_RT = gibbs_RT(); - for (int k = 0; k < m_kk; k++) { - xx = fmaxx(SmallNumber, m_s->moleFraction(k)); - mu[k] = rt*(g_RT[k] + log(xx) + logp); - } - } - - // new methods defined here - - - void IdealGasThermo::initThermo(Phase& s) { - Thermo::initThermo(s); - m_kk = s.nSpecies(); - m_mm = s.nElements(); - doublereal tmin = m_spthermo->minTemp(); - doublereal tmax = m_spthermo->maxTemp(); - if (tmin > 0.0) m_tmin = tmin; - if (tmax > 0.0) m_tmax = tmax; - m_p0 = refPressure(); - - // allocate space to cache species thermo properties - m_kk = m_s->nSpecies(); - - int leng = m_kk; - m_h0_RT.resize(leng); - m_g0_RT.resize(leng); - m_expg0_RT.resize(leng); - m_cp0_R.resize(leng); - m_s0_R.resize(leng); - m_pe.resize(leng, 0.0); - m_pp.resize(leng); - } - - - /** - * Set mixture to an equilibrium state consistent with specified - * element potentials and temperature. - * - * @param lambda_RT vector of non-dimensional element potentials - * \f[ \lambda_m/RT \f]. - * @param t temperature in K. - * @param work. Temporary work space. Must be dimensioned at least - * as large as the number of species. - * - */ - void IdealGasThermo::setToEquilState(const doublereal* lambda_RT) - { - const array_fp& grt = gibbs_RT(); - - // set the pressure and composition to be consistent with - // the temperature, - doublereal pres = 0.0; - for (int k = 0; k < m_kk; k++) { - m_pp[k] = -grt[k]; - for (int m = 0; m < m_mm; m++) { - m_pp[k] += phase().nAtoms(k,m)*lambda_RT[m]; - } - m_pp[k] = m_p0 * exp(m_pp[k]); - pres += m_pp[k]; - } - // set state - setState_PX(pres, m_pp.begin()); - } - - void IdealGasThermo::_updateThermo() const { - doublereal tnow = m_s->temperature(); - if (m_tlast != tnow) { - m_spthermo->update(tnow, m_cp0_R.begin(), m_h0_RT.begin(), - m_s0_R.begin()); - m_tlast = tnow; - doublereal rrt = 1.0 / (GasConstant * tnow); - int k; - doublereal deltaE; - for (k = 0; k < m_kk; k++) { - deltaE = rrt * m_pe[k]; - m_h0_RT[k] += deltaE; - m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; - } - m_tlast = tnow; - } - } -} - - - - diff --git a/Cantera/src/ImplicitChem.cpp b/Cantera/src/ImplicitChem.cpp deleted file mode 100755 index ea2760aa0..000000000 --- a/Cantera/src/ImplicitChem.cpp +++ /dev/null @@ -1,88 +0,0 @@ -/** - * @file Reactor.cpp - */ - -/* $Author$ - * $Revision$ - * $Date$ - */ - -// Copyright 2001 California Institute of Technology - - -#ifdef WIN32 -#pragma warning(disable:4786) -#pragma warning(disable:4503) -#endif - -#include "ImplicitChem.h" -#include "CVode.h" - -namespace Cantera { - - ImplicitChem::ImplicitChem(Kinetics& kin, ThermoPhase& therm) - : FuncEval(), m_kin(&kin), m_thermo(&therm), m_integ(0), - m_atol(1.e-15), m_rtol(1.e-7), m_maxstep(0.0), m_energy(false) - { - m_integ = new CVodeInt; - //m_mix = &kin.phase(); - m_wt = m_thermo->molecularWeights(); - - // use backward differencing, with a full Jacobian computed - // numerically, and use a Newton linear iterator - m_integ->setMethod(BDF_Method); - m_integ->setProblemType(DENSE + NOJAC); - m_integ->setIterator(Newton_Iter); - m_nsp = m_thermo->nSpecies(); - } - - // overloaded method of FuncEval. Called by the integrator to - // get the initial conditions. - void ImplicitChem::getInitialConditions(double t0, size_t leny, double* y) - { - m_thermo->getMassFractions(y); - m_h0 = m_thermo->enthalpy_mass(); - m_rho = m_thermo->density(); - m_press = m_thermo->pressure(); - } - - - /** - * Must be called before calling method 'advance' - */ - void ImplicitChem::initialize(doublereal t0) { - m_integ->setTolerances(m_rtol, m_atol); - // m_integ->setMaxStep(m_maxstep); - m_integ->initialize(t0, *this); - } - - - void ImplicitChem::updateState(doublereal* y) { - m_thermo->setMassFractions(y); - if (m_energy) { - doublereal delta, temp = m_thermo->temperature(); - do { - delta = -(m_thermo->enthalpy_mass() - m_h0)/m_thermo->cp_mass(); - temp += delta; - m_thermo->setTemperature(temp); - } - while (fabs(delta) > 1.e-7); - } - m_thermo->setPressure(m_press); - } - - /** - * Called by the integrator to evaluate ydot given y at time 'time'. - */ - void ImplicitChem::eval(doublereal time, doublereal* y, doublereal* ydot) - { - updateState(y); // synchronize the mixture state with y - m_thermo->setPressure(m_press); - m_kin->getNetProductionRates(ydot); // "omega dot" - int k; - for (k = 0; k < m_nsp; k++) { - ydot[k] *= m_wt[k]/m_rho; - } - } - -} diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in index 7d26212a1..f28e18ee0 100755 --- a/Cantera/src/Makefile.in +++ b/Cantera/src/Makefile.in @@ -54,7 +54,7 @@ FLOW1D = $(KINETICS) $(SOLVERS) EVERYTHING = $(KINETICS) $(HETEROKIN) $(ELECTROCHEM) $(EQUIL) $(CK) \ $(TRANSPORT) $(REACTOR) $(RPATH) $(SOLVERS) $(FLOW1D) -PCH = all.h +PCH = #config.h ct_defs.h utilities.h ThermoPhase.h Kinetics.h ReactionData.h RateCoeffMgr.h ReactionStoichMgr.h PCHGCH = $(PCH:.h=.h.gch) diff --git a/Cantera/src/SolidCompound.cpp b/Cantera/src/SolidCompound.cpp deleted file mode 100644 index c7f7ca9dd..000000000 --- a/Cantera/src/SolidCompound.cpp +++ /dev/null @@ -1,50 +0,0 @@ -/** - * - * @file StoichSubstance.cpp - * - */ - -#ifdef WIN32 -#pragma warning(disable:4786) -#pragma warning(disable:4503) -#endif - -#include "ct_defs.h" -#include "mix_defs.h" -#include "StoichSubstance.h" -#include "SpeciesThermo.h" - -namespace Cantera { - - void StoichSubstance::initThermo() { - m_kk = nSpecies(); - if (m_kk > 1) { - throw CanteraError("initThermo", - "stoichiometric substances may only contain one species."); - } - doublereal tmin = m_spthermo->minTemp(); - doublereal tmax = m_spthermo->maxTemp(); - if (tmin > 0.0) m_tmin = tmin; - if (tmax > 0.0) m_tmax = tmax; - m_p0 = refPressure(); - - int leng = m_kk; - m_h0_RT.resize(leng); - m_cp0_R.resize(leng); - m_s0_R.resize(leng); - } - - - void StoichSubstance::_updateThermo() const { - doublereal tnow = temperature(); - if (m_tlast != tnow) { - m_spthermo->update(tnow, m_cp0_R.begin(), m_h0_RT.begin(), - m_s0_R.begin()); - m_tlast = tnow; - } - } -} - - - - diff --git a/Cantera/src/importSurfChem.cpp b/Cantera/src/importSurfChem.cpp deleted file mode 100755 index b2be5b8a4..000000000 --- a/Cantera/src/importSurfChem.cpp +++ /dev/null @@ -1,106 +0,0 @@ - -#include "surfKinetics.h" -#include "ctml.h" -using namespace ctml; - -namespace Cantera { - - - /** - * Import a surface reaction mechanism - */ - void importInterfaceData(SurfacePhase* ph, SurfKinetics* kin, - string fname, string id) { - - ifstream f(fname.c_str()); - XML_Node root; - root.build(f); - - XML_Node* srxns = root.findID(id); - map fmap; - getFloats(*srxns, fmap); - ph->setSiteDensity(fmap["site_density"]); - XML_Node& spset = srxns->child("SpeciesArray"); - vector sp; - spset.getChildren("species",sp); - int nsp = sp.size(); - int k; - for (k = 0; k < nsp; k++) { - XML_Node& s = *sp[k]; - ph->addSpecies(s["name"], atof(s["size"].c_str())); - } - - vector rxns; - srxns->child("ReactionArray").getChildren("reaction",rxns); - int nrxns = rxns.size(); - int i, n; - string phase; - vector_int rindex, order, rstoich, pindex, pstoich; - - // get bulk phase data - int kk1 = kin->bulkPhase(0)->nSpecies(); - int kk2 = 0; - if (kin->bulkPhase(1)) kk2 = kin->bulkPhase(1)->nSpecies(); - vector bphase; - srxns->getChildren("phase", bphase); - int nbulk = bphase.size(); - string s, t; - vector phase_id(2,""); - for (int nb = 0; nb < nbulk; nb++) { - phase_id[nb] = (*bphase[nb])["id"]; - } - for (i = 0; i < nrxns; i++) { - XML_Node& rxn = *rxns[i]; - vector reac; - rxn.getChildren("reactant",reac); - int nr = reac.size(); - int k; - for (n = 0; n < nr; n++) { - XML_Node& r = *reac[n]; - rstoich.push_back(atoi(r["stoich"].c_str())); - order.push_back(atoi(r["order"].c_str())); - phase = r["phase"]; - if (phase == phase_id[0]) { - k = kin->bulkPhase(0)->speciesIndex(r["name"]); - } - else if (phase == phase_id[1]) { - k = kin->bulkPhase(1)->speciesIndex(r["name"]) + kk1; - } - else { - k = ph->speciesIndex(r["name"]) + kk1 + kk2; - } - rindex.push_back(k); - } - - vector prod; - rxn.getChildren("product",prod); - int np = prod.size(); - for (n = 0; n < np; n++) { - XML_Node& p = *prod[n]; - pstoich.push_back(atoi(p["stoich"].c_str())); - phase = p["phase"]; - if (phase == phase_id[0]) { - k = kin->bulkPhase(0)->speciesIndex(p["name"]); - } - else if (phase == phase_id[1]) { - k = kin->bulkPhase(1)->speciesIndex(p["name"]) + kk1; - } - else { - k = ph->speciesIndex(p["name"]) + kk1 + kk2; - } - pindex.push_back(k); - } - - XML_Node& rate = rxn.child("rate"); - map rp; - getFloats(rate, rp); - vector_fp kf(3); - kf[0] = rp["A"]; - kf[1] = rp["n"]; - kf[2] = rp["E"]; - - kin->addReaction(rindex, rstoich, order, pindex, pstoich, kf); - } - } - -} diff --git a/Cantera/src/surfKinetics.cpp b/Cantera/src/surfKinetics.cpp deleted file mode 100755 index afd72f7d7..000000000 --- a/Cantera/src/surfKinetics.cpp +++ /dev/null @@ -1,391 +0,0 @@ -/** - * @file SurfKinetics.cpp - * - */ - -// Copyright 2002 California Institute of Technology - - -// turn off warnings under Windows -#ifdef WIN32 -#pragma warning(disable:4786) -#pragma warning(disable:4503) -#endif - - -#include "surfKinetics.h" -#include "ReactionData.h" -#include "RateCoeffMgr.h" -#include "ImplicitSurfChem.h" -#include -using namespace std; - -#include "ctml.h" -using namespace ctml; - -#include - -namespace Cantera { - - void importInterfaceData(SurfacePhase* ph, SurfKinetics* kin, string fname, string id); - - /** - * Construct an empty surface reaction mechanism. - */ - SurfKinetics:: - SurfKinetics(SurfacePhase* surfphase, - thermo_t* th1, - thermo_t* th2, string fname, string id) : - Kinetics(), - m_surfphase(surfphase), - m_kk(0), - m_kk1(0), - m_kk2(0), - m_ktot(0), - m_nirrev(0), - m_integrator(0), - m_finalized(false), - m_twobulk(false), - m_xml(new XML_Node("interface_reactions")) - { - - // add the two bulk phases - addPhase(*th1); - if (th2) { - m_twobulk = true; - addPhase(*th2); - } - - m_kk1 = phase(0).nSpecies(); - - if (th2) { - m_kk2 = phase(1).nSpecies(); - } - m_kk = m_surfphase->nSpecies(); - m_kdata = new SurfKineticsData; - m_kdata->m_temp = 0.0; - - if (fname != "") importInterfaceData(surfphase, this, fname, id); - } - - void SurfKinetics:: - _update_rates_T() { - doublereal T = m_surfphase->temperature(); - if (T != m_kdata->m_temp) { - doublereal logT = log(T); - m_rates.update(T, logT, m_kdata->m_rfn.begin()); - m_kdata->m_temp = T; - m_kdata->m_ROP_ok = false; - } - }; - - void SurfKinetics:: - _update_rates_C() { - phase(0).getConcentrations(m_conc.begin()); - if (m_twobulk) { - phase(1).getConcentrations(m_conc.begin() + m_kk1); - } - m_surfphase->getConcentrations(m_conc.begin() + m_kk1 + m_kk2); - m_rates.update_C(m_conc.begin()); - m_kdata->m_ROP_ok = false; - } - - void SurfKinetics::updateROP() { - - _update_rates_C(); - _update_rates_T(); - - if (m_kdata->m_ROP_ok) return; - - const vector_fp& rf = m_kdata->m_rfn; - vector_fp& ropf = m_kdata->m_ropf; - - // copy rate coefficients into ropf - copy(rf.begin(), rf.end(), ropf.begin()); - - // multiply by perturbation factor - multiply_each(ropf.begin(), ropf.end(), m_perturb.begin()); - - // multiply ropf by concentration products - int i, j, k, o; - for (i = 0; i < m_ii; i++) { - for (j = 0; j < m_nr[i]; j++) { - k = m_reactants[i][j]; - o = m_order[i][j]; - ropf[i] *= pow(m_conc[k],m_order[i][j]); - } - } - m_kdata->m_ROP_ok = true; - } - - - void SurfKinetics:: - getNetProductionRates(doublereal* net) { - updateROP(); - int i, n, k; - doublereal q; - for (k = 0; k < m_ktot; k++) net[k] = 0.0; - - for (i = 0; i < m_ii; i++) { - q = m_kdata->m_ropf[i]; - for (n = 0; n < m_nr[i]; n++) { - k = m_reactants[i][n]; - net[k] -= q*m_rst[i][n]; - } - for (n = 0; n < m_np[i]; n++) { - k = m_products[i][n]; - net[k] += q*m_pst[i][n]; - } - } - } - - void SurfKinetics:: - getCreationRates(doublereal* cdot) { - updateROP(); - int i, n, k; - doublereal q; - fill(cdot, cdot + m_ktot, 0.0); - for (i = 0; i < m_ii; i++) { - q = m_kdata->m_ropf[i]; - for (n = 0; n < m_np[i]; n++) { - k = m_products[i][n]; - cdot[k] += q*m_pst[i][n]; - } - } - } - - void SurfKinetics:: - getDestructionRates(doublereal* ddot) { - updateROP(); - int i, n, k; - doublereal q; - fill(ddot, ddot + m_ktot, 0.0); - for (i = 0; i < m_ii; i++) { - q = m_kdata->m_ropf[i]; - for (n = 0; n < m_nr[i]; n++) { - k = m_reactants[i][n]; - ddot[k] += q*m_rst[i][n]; - } - } - } - - void SurfKinetics:: - getChemRates(doublereal* rtau) { - updateROP(); - int i, n, k; - doublereal q; - fill(rtau, rtau + m_ktot, 0.0); - for (i = 0; i < m_ii; i++) { - q = m_kdata->m_ropf[i]; - for (n = 0; n < m_nr[i]; n++) { - k = m_reactants[i][n]; - rtau[k] += q*m_rst[i][n]; - } - } - for (k = 0;k < m_ktot; k++) { - if (m_conc[k] != 0.0) - rtau[k] = fabs(rtau[k]/m_conc[k]); - else - rtau[k] = 0.0; - } - } - - void SurfKinetics:: - saveReactionData( - const vector_int& r, - const vector_int& rstoich, - const vector_int& order, - const vector_int& p, - const vector_int& pstoich, - const vector_fp& rateParams) { - - if (nReactions() == 0) - m_xml->addChild("ReactionArray"); - - XML_Node& rxndata = *new XML_Node("reaction"); - int n, k; - string nm, ph, ustr, comment; - for (n = 0; n < r.size(); n++) { - XML_Node& reac = rxndata.addChild("reactant"); - if (r[n] < m_kk1) { - k = r[n]; - nm = phase(0).speciesName(k); - ph = phase(0).id(); - ustr = "kmol/m^3"; - m_bsp1[nm] = 1; - } - else if (r[n] < m_kk1 + m_kk2) { - k = r[n] - m_kk1; - nm = phase(1).speciesName(k); - ph = phase(1).id(); - ustr = "kmol/m^3"; - m_bsp2[nm] = 1; - } - else { - k = r[n] - m_kk1 - m_kk2; nm = m_surfphase->speciesName(k); - ph = ""; // m_surfphase->id(); - ustr = "kmol/m^2"; - } - if (ph != "") reac.addAttribute("phase",ph); - reac.addAttribute("name",nm); - reac.addAttribute("stoich",rstoich[n]); - reac.addAttribute("order",order[n]); - // reac.addAttribute("units",ustr); - comment += nm+" + "; - } - comment = comment.substr(0, comment.size() - 2) + " => "; - - for (n = 0; n < p.size(); n++) { - XML_Node& prod = rxndata.addChild("product"); - if (p[n] < m_kk1) { - k = p[n]; nm = phase(0).speciesName(k); - ph = phase(0).id(); - ustr = "kmol/m^3"; - } - else if (p[n] < m_kk1 + m_kk2) { - k = p[n] - m_kk1; - nm = phase(1).speciesName(k); - ph = phase(1).id(); - ustr = "kmol/m^3"; - } - else { - k = p[n] - m_kk1 - m_kk2; - nm = m_surfphase->speciesName(k); - ph = ""; - ustr = "kmol/m^2"; - } - if (ph != "") prod.addAttribute("phase",ph); - prod.addAttribute("name",nm); - prod.addAttribute("stoich",pstoich[n]); - comment += nm+" + "; - } - comment = " "+comment.substr(0, comment.size() - 2)+" "; - - XML_Node& rate = rxndata.addChild("rate"); - rate.addAttribute("type","Arrhenius"); - rate.addAttribute("units","kmol/m^2/s"); - addFloat(rate, "A", rateParams[0]); - addFloat(rate, "n", rateParams[1]); - addFloat(rate, "E", rateParams[2], "K"); - - XML_Node& rxns = m_xml->child("ReactionArray"); - rxns.addComment(comment); - rxns.addChild(rxndata); - }; - - - void SurfKinetics:: - addReaction(const vector_int& r, - const vector_int& rstoich, - const vector_int& order, - const vector_int& p, - const vector_int& pstoich, - const vector_fp& rateParams) { - - // record reaction parameters - saveReactionData(r, rstoich, order, p, pstoich, rateParams); - - // prohibit adding more species - if (!m_surfphase->speciesFrozen()) - m_surfphase->freezeSpecies(); - - // if init() hasn't been called yet, call it - if (m_kk == 0) init(); - - int iloc; - // install rate coeff calculator - iloc = m_rates.install( m_ii, - ARRHENIUS, rateParams.size(), rateParams.begin()); - - // add constant term to rate coeff value vector - m_kdata->m_rfn.push_back(rateParams[0]); - - // forward rxn order - m_order.push_back(order); - - m_kdata->m_ropf.push_back(0.0); // extend by one for new rxn - - m_reactants.push_back(r); - m_rst.push_back(rstoich); - m_products.push_back(p); - m_pst.push_back(pstoich); - - m_nr.push_back(r.size()); - m_np.push_back(p.size()); - - incrementRxnCount(); - } - - - void SurfKinetics::init() { - m_kk = m_surfphase->nSpecies(); - m_ktot = m_kk + m_kk1 + m_kk2; - m_conc.resize(m_ktot); - Kinetics::init(); - } - - void SurfKinetics::save(string fname, string idtag, string comment) { - struct tm *newtime; - time_t aclock; - ::time( &aclock ); /* Get time in seconds */ - newtime = localtime( &aclock ); /* Convert time to struct tm form */ - - ofstream fout(fname.c_str()); - XML_Node root("doc"); - XML_Node& ct = root.addChild("ctml"); - ct.addComment(comment); - - XML_Node& iface = ct.addChild("interface"); - addString(iface,"timestamp",asctime(newtime)); - iface.addAttribute("id",idtag); - addFloat(iface, "site_density", m_surfphase->siteDensity()); - XML_Node& bp1 = iface.addChild("phase"); - bp1.addAttribute("id",phase(0).id()); - map::const_iterator b = m_bsp1.begin(), e = m_bsp1.end(); - for (; b != e; ++b) { - bp1.addChild("species").addAttribute("name",b->first); - } - bp1.addChild(thermo(0).xml()); - if (m_twobulk) { - XML_Node& bp2 = iface.addChild("phase"); - bp2.addAttribute("id",phase(1).id()); - map::const_iterator b = m_bsp2.begin(), e = m_bsp2.end(); - for (; b != e; ++b) { - bp2.addChild("species").addAttribute("name",b->first); - } - bp2.addChild(thermo(1).xml()); - } - iface.addChild(m_surfphase->xml().child("SpeciesArray")); - iface.addChild(m_xml->child("ReactionArray")); - ct.writeHeader(fout); - ct.write(fout); - fout.close(); - } - - void SurfKinetics::finalize() { - if (!m_finalized) { - m_finalized = true; - } - } - - bool SurfKinetics::ready() const { - return (m_finalized); - } - - void SurfKinetics::integrate(doublereal dt) { - finalize(); - if (m_integrator == 0) { - m_integrator = new ImplicitSurfChem(*this); - m_integrator->initialize(0.0); - } - m_integrator->integrate(0.0, dt); - } -} - - - - - - - -