cleanup
This commit is contained in:
Dave Goodwin
parent
3e7e28b6fb
commit
1ae73ae1c5
6 changed files with 1 additions and 744 deletions
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -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<string, double> fmap;
|
||||
getFloats(*srxns, fmap);
|
||||
ph->setSiteDensity(fmap["site_density"]);
|
||||
XML_Node& spset = srxns->child("SpeciesArray");
|
||||
vector<XML_Node*> 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<XML_Node*> 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<XML_Node*> bphase;
|
||||
srxns->getChildren("phase", bphase);
|
||||
int nbulk = bphase.size();
|
||||
string s, t;
|
||||
vector<string> phase_id(2,"<none>");
|
||||
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<XML_Node*> 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<XML_Node*> 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<string, doublereal> 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);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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 <iostream>
|
||||
using namespace std;
|
||||
|
||||
#include "ctml.h"
|
||||
using namespace ctml;
|
||||
|
||||
#include <time.h>
|
||||
|
||||
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<string,int>::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<string,int>::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);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
Add table
Reference in a new issue