research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Added a capability to solve the surface pseudosteady state problem

Browse source 
via a damped newton's method. It also employs a false transient
algorithm when far from the solution.
This commit is contained in:
Harry Moffat 2007-08-23 21:43:07 +00:00
parent 12e9694043
commit 60a93559af
7 changed files with 3115 additions and 531 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

368
Cantera/src/kinetics/ImplicitSurfChem.cpp
View file

@ -21,72 +21,135 @@
#include "ImplicitSurfChem.h"
#include "Integrator.h"
#include "solveSP.h"
using namespace std;
namespace Cantera {
// Constructor
ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
FuncEval(),
m_nsurf(0),
m_nv(0),
m_numBulkPhases(0),
m_numTotalBulkSpecies(0),
m_numTotalSpecies(0),
m_integ(0),
m_atol(1.e-14),
m_rtol(1.e-7),
m_maxstep(0.0),
m_mediumSpeciesStart(-1),
m_bulkSpeciesStart(-1),
m_surfSpeciesStart(-1),
m_surfSolver(0),
m_commonTempPressForPhases(true),
m_ioFlag(0)
{
m_nsurf = static_cast<int>(k.size());
int ns, nsp;
int nt, ntmax = 0;
int kinSpIndex = 0;
// Loop over the number of surface kinetics objects
for (int n = 0; n < m_nsurf; n++) {
InterfaceKinetics *kinPtr = k[n];
m_vecKinPtrs.push_back(kinPtr);
ns = k[n]->surfacePhaseIndex();
if (ns < 0)
throw CanteraError("ImplicitSurfChem",
"kinetics manager contains no surface phase");
m_surfindex.push_back(ns);
m_surf.push_back((SurfPhase*)&k[n]->thermo(ns));
nsp = m_surf.back()->nSpecies();
m_nsp.push_back(nsp);
m_nv += m_nsp.back();
nt = k[n]->nTotalSpecies();
if (nt > ntmax) ntmax = nt;
m_specStartIndex.push_back(kinSpIndex);
kinSpIndex += nsp;
ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k)
: FuncEval(), m_nv(0), m_integ(0),
m_atol(1.e-14), m_rtol(1.e-7), m_maxstep(0.0)
{
m_nsurf = static_cast<int>(k.size());
int ns;
int nt, ntmax = 0;
for (int n = 0; n < m_nsurf; n++) {
m_kin.push_back(k[n]);
ns = k[n]->surfacePhaseIndex();
if (ns < 0)
throw CanteraError("ImplicitSurfChem",
"kinetics manager contains no surface phase");
m_surfindex.push_back(ns);
m_surf.push_back((SurfPhase*)&k[n]->thermo(ns));
m_nsp.push_back(m_surf.back()->nSpecies());
m_nv += m_nsp.back();
nt = k[n]->nTotalSpecies();
if (nt > ntmax) ntmax = nt;
}
m_integ = newIntegrator("CVODE");// CVodeInt;
// 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_work.resize(ntmax);
int nPhases = kinPtr->nPhases();
vector_int pLocTmp(nPhases);
int imatch = -1;
for (int ip = 0; ip < nPhases; ip++) {
if (ip != ns) {
ThermoPhase *thPtr = & kinPtr->thermo(ip);
if ((imatch = checkMatch(m_bulkPhases, thPtr)) < 0) {
m_bulkPhases.push_back(thPtr);
m_numBulkPhases++;
nsp = thPtr->nSpecies();
m_nspBulkPhases.push_back(nsp);
m_numTotalBulkSpecies += nsp;
imatch = m_bulkPhases.size() - 1;
}
pLocTmp[ip] = imatch;
} else {
pLocTmp[ip] = -n;
}
}
pLocVec.push_back(pLocTmp);
}
m_numTotalSpecies = m_nv + m_numTotalBulkSpecies;
m_concSpecies.resize(m_numTotalSpecies, 0.0);
m_concSpeciesSave.resize(m_numTotalSpecies, 0.0);
m_integ = newIntegrator("CVODE");
/**
* Destructor. Deletes the integrator.
*/
ImplicitSurfChem::~ImplicitSurfChem(){
delete m_integ;
// 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_work.resize(ntmax);
}
// overloaded method of FuncEval. Called by the integrator to
// get the initial conditions.
void ImplicitSurfChem::getInitialConditions(double t0, size_t lenc,
double* c)
{
int loc = 0;
for (int n = 0; n < m_nsurf; n++) {
m_surf[n]->getCoverages(c + loc);
loc += m_nsp[n];
}
int ImplicitSurfChem::checkMatch(std::vector<ThermoPhase *> m_vec, ThermoPhase *thPtr) {
int retn = -1;
for (int i = 0; i < (int) m_vec.size(); i++) {
ThermoPhase *th = m_vec[i];
if (th == thPtr) {
return i;
}
}
return retn;
}
/*
* Destructor. Deletes the integrator.
*/
ImplicitSurfChem::~ImplicitSurfChem(){
if (m_integ) {
delete m_integ;
}
if (m_surfSolver) {
delete m_surfSolver;
}
}
// overloaded method of FuncEval. Called by the integrator to
// get the initial conditions.
void ImplicitSurfChem::getInitialConditions(double t0, size_t lenc,
double* c)
{
int loc = 0;
for (int n = 0; n < m_nsurf; n++) {
m_surf[n]->getCoverages(c + loc);
loc += m_nsp[n];
}
}
/**
* Must be called before calling method 'advance'
*/
void ImplicitSurfChem::initialize(doublereal t0) {
m_integ->setTolerances(m_rtol, m_atol);
m_integ->initialize(t0, *this);
}
/*
* Must be called before calling method 'advance'
*/
void ImplicitSurfChem::initialize(doublereal t0) {
m_integ->setTolerances(m_rtol, m_atol);
m_integ->initialize(t0, *this);
}
// Integrate from t0 to t1. The integrator is reinitialized first.
/*
@ -117,16 +180,15 @@ namespace Cantera {
updateState(m_integ->solution());
}
void ImplicitSurfChem::updateState(doublereal* c) {
int loc = 0;
for (int n = 0; n < m_nsurf; n++) {
m_surf[n]->setCoverages(c + loc);
loc += m_nsp[n];
}
void ImplicitSurfChem::updateState(doublereal* c) {
int loc = 0;
for (int n = 0; n < m_nsurf; n++) {
m_surf[n]->setCoverages(c + loc);
loc += m_nsp[n];
}
}
/**
/*
* Called by the integrator to evaluate ydot given y at time 'time'.
*/
void ImplicitSurfChem::eval(doublereal time, doublereal* y,
@ -138,8 +200,8 @@ namespace Cantera {
int loc, k, kstart;
for (n = 0; n < m_nsurf; n++) {
rs0 = 1.0/m_surf[n]->siteDensity();
m_kin[n]->getNetProductionRates(DATA_PTR(m_work));
kstart = m_kin[n]->kineticsSpeciesIndex(0,m_surfindex[n]);
m_vecKinPtrs[n]->getNetProductionRates(DATA_PTR(m_work));
kstart = m_vecKinPtrs[n]->kineticsSpeciesIndex(0,m_surfindex[n]);
sum = 0.0;
loc = 0;
for (k = 1; k < m_nsp[n]; k++) {
@ -151,4 +213,182 @@ namespace Cantera {
}
}
// Solve for the pseudo steady-state of the surface problem
/*
* Solve for the steady state of the surface problem.
* This is the same thing as the advanceCoverages() function,
* but at infinite times.
*
* Note, a direct solve is carried out under the hood here,
* to reduce the computational time.
*/
void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
doublereal timeScaleOverride) {
int ifunc;
/*
* set bulkFunc
* -> We assume that the bulk concentrations are constant.
*/
int bulkFunc = BULK_ETCH;
/*
* time scale - time over which to integrate equations
*/
double time_scale = timeScaleOverride;
/*
*
*/
if (!m_surfSolver) {
m_surfSolver = new solveSP(this, bulkFunc);
/*
* set ifunc, which sets the algorithm.
*/
ifunc = SFLUX_INITIALIZE;
} else {
ifunc = SFLUX_RESIDUAL;
}
// Possibly override the ifunc value
if (ifuncOverride >= 0) {
ifunc = ifuncOverride;
}
/*
* Get the specifications for the problem from the values
* in the ThermoPhase objects for all phases.
*
* 1) concentrations of all species in all phases, m_concSpecies[]
* 2) Temperature and pressure
*/
getConcSpecies(DATA_PTR(m_concSpecies));
InterfaceKinetics *ik = m_vecKinPtrs[0];
ThermoPhase &tp = ik->thermo(0);
double TKelvin = tp.temperature();
double PGas = tp.pressure();
/*
* Make sure that there is a common temperature and
* pressure for all ThermoPhase objects belonging to the
* interfacial kinetics object, if it is required by
* the problem statement.
*/
if (m_commonTempPressForPhases) {
setCommonState_TP(TKelvin, PGas);
}
double reltol = 1.0E-6;
double atol = 1.0E-20;
/*
* Install a filter for negative concentrations. One of the
* few ways solvess can fail is if concentrations on input
* are below zero.
*/
bool rset = false;
for (int k = 0; k < m_nv; k++) {
if (m_concSpecies[k] < 0.0) {
rset = true;
m_concSpecies[k] = 0.0;
}
}
if (rset) {
setConcSpecies(DATA_PTR(m_concSpecies));
}
m_surfSolver->ioflag = m_ioFlag;
// Save the current solution
copy(m_concSpecies.begin(), m_concSpecies.end(), m_concSpeciesSave.begin());
int retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
reltol, atol);
if (retn != 1) {
// reset the concentrations
copy(m_concSpeciesSave.begin(), m_concSpeciesSave.end(), m_concSpecies.begin());
setConcSpecies(DATA_PTR(m_concSpeciesSave));
ifunc = SFLUX_INITIALIZE;
retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
reltol, atol);
if (retn != 1) {
throw CanteraError("ImplicitSurfChem::solvePseudoSteadyStateProblem",
"solveSP return an error condition!");
}
}
}
/*
* getConcSpecies():
*
* Fills the local concentration vector, m_concSpecies for all of the
* species in all of the phases that are unknowns in the surface
* problem.
*
* m_concSpecies[]
*/
void ImplicitSurfChem::getConcSpecies(doublereal * const vecConcSpecies) const {
int kstart;
for (int ip = 0; ip < m_nsurf; ip++) {
ThermoPhase * TP_ptr = m_surf[ip];
kstart = m_specStartIndex[ip];
TP_ptr->getConcentrations(vecConcSpecies + kstart);
}
kstart = m_nv;
for (int ip = 0; ip < m_numBulkPhases; ip++) {
ThermoPhase * TP_ptr = m_bulkPhases[ip];
int nsp = TP_ptr->nSpecies();
TP_ptr->getConcentrations(vecConcSpecies + kstart);
kstart += nsp;
}
}
/*
* setConcSpecies():
*
* Fills the local concentration vector, m_concSpecies for all of the
* species in all of the phases that are unknowns in the surface
* problem.
*
* m_concSpecies[]
*/
void ImplicitSurfChem::setConcSpecies(const doublereal * const vecConcSpecies) {
int kstart;
for (int ip = 0; ip < m_nsurf; ip++) {
ThermoPhase * TP_ptr = m_surf[ip];
kstart = m_specStartIndex[ip];
TP_ptr->setConcentrations(vecConcSpecies + kstart);
}
kstart = m_nv;
for (int ip = 0; ip < m_numBulkPhases; ip++) {
ThermoPhase * TP_ptr = m_bulkPhases[ip];
int nsp = TP_ptr->nSpecies();
TP_ptr->setConcentrations(vecConcSpecies + kstart);
kstart += nsp;
}
}
/*
* setCommonState_TP():
*
* Sets a common temperature and pressure amongst the
* thermodynamic objects in the interfacial kinetics object.
*
* Units Temperature = Kelvin
* Pressure = Pascal
*/
void ImplicitSurfChem::
setCommonState_TP(double TKelvin, double PresPa) {
int nphases = m_nsurf;
for (int ip = 0; ip < nphases; ip++) {
ThermoPhase *TP_ptr = m_surf[ip];
TP_ptr->setState_TP(TKelvin, PresPa);
}
for (int ip = 0; ip < m_numBulkPhases; ip++) {
ThermoPhase *TP_ptr = m_bulkPhases[ip];
TP_ptr->setState_TP(TKelvin, PresPa);
}
}
}

176
Cantera/src/kinetics/ImplicitSurfChem.h
View file

@ -26,9 +26,12 @@
#include "Integrator.h"
#include "InterfaceKinetics.h"
#include "SurfPhase.h"
#include "solveSP.h"
namespace Cantera {
class solveSP;
//! Advances the surface coverages of the associated set of SurfacePhase
//! objects in time
@ -115,6 +118,36 @@ namespace Cantera {
* @param t1 Final Time -> This is an input
*/
void integrate0(doublereal t0, doublereal t1);
//! Solve for the pseudo steady-state of the surface problem
/*!
* Solve for the steady state of the surface problem.
* This is the same thing as the advanceCoverages() function,
* but at infinite times.
*
* Note, a direct solve is carried out under the hood here,
* to reduce the computational time.
*
* @param ifuncOverride 4 values are possible
* 1 SFLUX_INITIALIZE
* 2 SFLUX_RESIDUAL
* 3 SFLUX_JACOBIAN
* 4 SFLUX_TRANSIENT
* The default is -1, which means that the program
* will decide.
* @param timeScaleOverride When a psuedo transient is
* selected this value can be used to override
* the default time scale for integration which
* is one.
* When SFLUX_TRANSIENT is used, this is equal to the
* time over which the equations are integrated.
* When SFLUX_INITIALIZE is used, this is equal to the
* time used in the initial transient algorithm,
* before the equation system is solved directly.
*/
void solvePseudoSteadyStateProblem(int ifuncOverride = -1,
doublereal timeScaleOverride = 1.0);
// overloaded methods of class FuncEval
@ -144,6 +177,59 @@ namespace Cantera {
virtual void getInitialConditions(doublereal t0,
size_t leny, doublereal* y);
/*
* Get the specifications for the problem from the values
* in the ThermoPhase objects for all phases.
*
* 1) concentrations of all species in all phases, m_concSpecies[]
* 2) Temperature and pressure
*
*
* @param vecConcSpecies Vector of concentrations. The
* phase concentration vectors are contiguous
* within the object, in the same order as the
* unknown vector.
*/
void getConcSpecies(doublereal * const vecConcSpecies) const;
//! Sets the concentrations within phases that are unknowns in
//! the surface problem
/*!
* Fills the local concentration vector for all of the
* species in all of the phases that are unknowns in the surface
* problem.
*
* @param vecConcSpecies Vector of concentrations. The
* phase concentration vectors are contiguous
* within the object, in the same order as the
* unknown vector.
*/
void setConcSpecies(const doublereal * const vecConcSpecies);
//! Sets the state variable in all thermodynamic phases (surface and
//! surrounding bulk phases) to the input temperature and pressure
/*!
* @param TKelvin input temperature (kelvin)
* @param PresPa input pressure in pascal.
*/
void setCommonState_TP(double TKelvin, double PresPa);
//! Returns a reference to the vector of pointers to the
//! InterfaceKinetics objects
/*!
* This should probably go away in the future, as it opens up the
* class.
*/
std::vector<InterfaceKinetics*> & getObjects() {
return m_vecKinPtrs;
}
int checkMatch(std::vector<ThermoPhase *> m_vec, ThermoPhase *thPtr);
void setIOFlag(int ioFlag) {
m_ioFlag = ioFlag;
}
protected:
@ -161,20 +247,100 @@ namespace Cantera {
*/
void updateState(doublereal* y);
//! vector of pointers to surface phases.
std::vector<SurfPhase*> m_surf;
std::vector<InterfaceKinetics*> m_kin;
//! Vector of pointers to bulk phases
std::vector<ThermoPhase *> m_bulkPhases;
//! vector of pointers to InterfaceKinetics objects
std::vector<InterfaceKinetics*> m_vecKinPtrs;
//! Vector of number of species in each Surface Phase
vector_int m_nsp;
//! index of the surface phase in each InterfaceKinetics object
vector_int m_surfindex;
int m_nsurf;
int m_nv;
//int m_nsp, m_surfindex;
Integrator* m_integ; // pointer to integrator
vector_int m_specStartIndex;
//! Total number of surface phases.
/*!
* This is also equal to the number of InterfaceKinetics objects
* as there is a 1-1 correspondence between InterfaceKinetics objects
* and surface phases.
*/
int m_nsurf;
//! Total number of surface species in all surface phases
/*!
* This is the total number of unknowns in m_mode 0 problem
*/
int m_nv;
int m_numBulkPhases;
vector_int m_nspBulkPhases;
int m_numTotalBulkSpecies;
int m_numTotalSpecies;
std::vector<vector_int> pLocVec;
//! Pointer to the cvode integrator
Integrator* m_integ;
doublereal m_atol, m_rtol; // tolerances
doublereal m_maxstep; // max step size
vector_fp m_work;
/**
* Temporary vector - length num species in the Kinetics object.
* This is the sum of the number of species
* in each phase included in the kinetics object.
*/
vector_fp m_concSpecies;
vector_fp m_concSpeciesSave;
//std::vector<vector_fp> m_vectorConcKinSpecies;
//std::vector<vector_fp> m_vectorNetSpeciesProdRate;
/**
* Index into the species vector of the kinetics manager,
* pointing to the first species from the surrounding medium.
*/
int m_mediumSpeciesStart;
/**
* Index into the species vector of the kinetics manager,
* pointing to the first species from the condensed phase
* of the particles.
*/
int m_bulkSpeciesStart;
/**
* Index into the species vector of the kinetics manager,
* pointing to the first species from the surface
* of the particles
*/
int m_surfSpeciesStart;
/**
* Pointer to the helper method, Placid, which solves the
* surface problem.
*/
solveSP *m_surfSolver;
//! If true, a common temperature and pressure for all
//! surface and bulk phases associated with the surface problem
//! is imposed
bool m_commonTempPressForPhases;
//! We make the solveSS class a friend because we need
//! to access all of the above information directly.
//! Adding the members into the class is also a possibility.
friend class solveSS;
private:
//! Controls the amount of printing from this routine
//! and underlying routines.
int m_ioFlag;
};
}

60
Cantera/src/kinetics/InterfaceKinetics.cpp
View file

@ -46,7 +46,8 @@ namespace Cantera {
m_integrator(0),
m_finalized(false),
m_has_coverage_dependence(false),
m_has_electrochem_rxns(false)
m_has_electrochem_rxns(false),
m_ioFlag(0)
{
if (thermo != 0) addPhase(*thermo);
m_kdata = new InterfaceKineticsData;
@ -59,7 +60,9 @@ namespace Cantera {
InterfaceKinetics::
~InterfaceKinetics(){
delete m_kdata;
delete m_integrator;
if (m_integrator) {
delete m_integrator;
}
}
@ -698,6 +701,13 @@ namespace Cantera {
}
void InterfaceKinetics::setIOFlag(int ioFlag) {
m_ioFlag = ioFlag;
if (m_integrator) {
m_integrator->setIOFlag(ioFlag);
}
}
// void InterfaceKinetics::
// addGlobalReaction(const ReactionData& r) {
@ -905,31 +915,39 @@ namespace Cantera {
*
* Note, a direct solve is carried out under the hood here,
* to reduce the computational time.
*
* the integrator object is saved inbetween calls to
* reduce the computational cost of repeated calls.
*/
void InterfaceKinetics::solvePseudoSteadyStateProblem() {
#ifndef DEBUG_HKM
advanceCoverages(1000.0);
#else
void InterfaceKinetics::
solvePseudoSteadyStateProblem(int ifuncOverride,
doublereal timeScaleOverride) {
// create our own solver object
if (m_integrator == 0) {
vector<InterfaceKinetics*> k;
k.push_back(this);
m_integrator = new ImplicitSurfChem(k);
m_integrator->initialize();
}
m_integrator->setIOFlag(m_ioFlag);
/*
* New direct method to go here
*/
#endif
m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
}
void EdgeKinetics::finalize() {
m_rwork.resize(nReactions());
int ks = reactionPhaseIndex();
if (ks < 0) throw CanteraError("EdgeKinetics::finalize",
"no edge phase is present.");
m_surf = (SurfPhase*)&thermo(ks);
if (m_surf->nDim() != 1)
throw CanteraError("EdgeKinetics::finalize",
"expected interface dimension = 1, but got dimension = "
+int2str(m_surf->nDim()));
m_finalized = true;
}
void EdgeKinetics::finalize() {
m_rwork.resize(nReactions());
int ks = reactionPhaseIndex();
if (ks < 0) throw CanteraError("EdgeKinetics::finalize",
"no edge phase is present.");
m_surf = (SurfPhase*)&thermo(ks);
if (m_surf->nDim() != 1)
throw CanteraError("EdgeKinetics::finalize",
"expected interface dimension = 1, but got dimension = "
+int2str(m_surf->nDim()));
m_finalized = true;
}
}

905
Cantera/src/kinetics/InterfaceKinetics.h
View file

@ -1,6 +1,8 @@
/**
* @file InterfaceKinetics.h
*
*/
/*
* $Author$
* $Revision$
* $Date$
@ -60,519 +62,546 @@ namespace Cantera {
};
///
/// A kinetics manager for heterogeneous reaction mechanisms. The
/// reactions are assumed to occur at a 2D interface between two
/// 3D phases.
///
class InterfaceKinetics : public Kinetics {
public:
/**
* Constructor
*
* @param thermo The optional parameter may be used to initialize
* the object with one ThermoPhase object.
* HKM Note -> Since the interface kinetics
* object will probably require multiple thermophase
* objects, this is probably not a good idea
* to have this parameter.
*/
InterfaceKinetics(thermo_t* thermo = 0);
/// Destructor.
virtual ~InterfaceKinetics();
virtual int ID() { return cInterfaceKinetics; }
virtual int type() { return cInterfaceKinetics; }
/**
* Set the electric potential in the nth phase
*
* @param n phase Index in this kinetics object.
* @param V Electric potential (volts)
*/
void setElectricPotential(int n, doublereal V) {
thermo(n).setElectricPotential(V);
m_redo_rates = true;
}
///
/// A kinetics manager for heterogeneous reaction mechanisms. The
/// reactions are assumed to occur at a 2D interface between two
/// 3D phases.
/// @name Reaction Rates Of Progress
///
class InterfaceKinetics : public Kinetics {
//@{
public:
//! Return the forward rates of progress for each reaction
/*!
* @param fwdROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getFwdRatesOfProgress(doublereal* fwdROP) {
updateROP();
std::copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP);
}
/**
* Constructor
*
* @param thermo The optional parameter may be used to initialize
* the object with one ThermoPhase object.
* HKM Note -> Since the interface kinetics
* object will probably require multiple thermophase
* objects, this is probably not a good idea
* to have this parameter.
*/
InterfaceKinetics(thermo_t* thermo = 0);
//! Return the reverse rates of progress for each reaction
/*!
* @param revROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getRevRatesOfProgress(doublereal* revROP) {
updateROP();
std::copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP);
}
//! Return the net rates of progress for each reaction
/*!
* @param netROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getNetRatesOfProgress(doublereal* netROP) {
updateROP();
std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP);
}
virtual void getEquilibriumConstants(doublereal* kc);
/// Destructor.
virtual ~InterfaceKinetics();
virtual void getDeltaGibbs( doublereal* deltaG);
virtual int ID() { return cInterfaceKinetics; }
virtual int type() { return cInterfaceKinetics; }
/**
* Set the electric potential in the nth phase
*
* @param n phase Index in this kinetics object.
* @param V Electric potential (volts)
*/
void setElectricPotential(int n, doublereal V) {
thermo(n).setElectricPotential(V);
m_redo_rates = true;
}
///
/// @name Reaction Rates Of Progress
///
//@{
//! Return the forward rates of progress for each reaction
/*!
* @param fwdROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getFwdRatesOfProgress(doublereal* fwdROP) {
updateROP();
std::copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP);
}
//! Return the reverse rates of progress for each reaction
/*!
* @param revROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getRevRatesOfProgress(doublereal* revROP) {
updateROP();
std::copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP);
}
//! Return the net rates of progress for each reaction
/*!
* @param netROP vector of rates of progress.
* length = number of reactions, Units are kmol m-2 s-1.
*/
virtual void getNetRatesOfProgress(doublereal* netROP) {
updateROP();
std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP);
}
virtual void getEquilibriumConstants(doublereal* kc);
virtual void getDeltaGibbs( doublereal* deltaG);
/**
* Return the vector of values for the reactions change in
* enthalpy.
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1
*/
virtual void getDeltaEnthalpy( doublereal* deltaH);
/**
* Return the vector of values for the reactions change in
* enthalpy.
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1
*/
virtual void getDeltaEnthalpy( doublereal* deltaH);
//! Return the vector of values for the change in
//! entropy due to each reaction
/*!
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1 Kelvin-1
*
* @param deltaS vector of Enthalpy changes
* Length = m_ii, number of reactions
*
*/
virtual void getDeltaEntropy(doublereal* deltaS);
//! Return the vector of values for the change in
//! entropy due to each reaction
/*!
* These values depend upon the concentration
* of the solution.
*
* units = J kmol-1 Kelvin-1
*
* @param deltaS vector of Enthalpy changes
* Length = m_ii, number of reactions
*
*/
virtual void getDeltaEntropy(doublereal* deltaS);
//! Return the vector of values for the reaction
//! standard state gibbs free energy change.
/*!
* These values don't depend upon the concentration
* of the solution.
*
* @param deltaG vector of rxn SS free energy changes
* units = J kmol-1
*/
virtual void getDeltaSSGibbs(doublereal* deltaG);
//! Return the vector of values for the reaction
//! standard state gibbs free energy change.
/*!
* These values don't depend upon the concentration
* of the solution.
*
* @param deltaG vector of rxn SS free energy changes
* units = J kmol-1
*/
virtual void getDeltaSSGibbs(doublereal* deltaG);
//! Return the vector of values for the change in the
//! standard state enthalpies of reaction.
/*!
* These values don't depend upon the concentration
* of the solution.
*
* @param deltaH vector of rxn SS enthalpy changes
* units = J kmol-1
*/
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
//! Return the vector of values for the change in the
//! standard state enthalpies of reaction.
/*!
* These values don't depend upon the concentration
* of the solution.
*
* @param deltaH vector of rxn SS enthalpy changes
* units = J kmol-1
*/
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
//! Return the vector of values for the change in the
//! standard state entropies for each reaction.
/*!
* These values don't depend upon the concentration
* of the solution.
*
* @param deltaS vector of rxn SS entropy changes
* units = J kmol-1 Kelvin-1
*/
virtual void getDeltaSSEntropy(doublereal* deltaS);
//! Return the vector of values for the change in the
//! standard state entropies for each reaction.
/*!
* These values don't depend upon the concentration
* of the solution.
*
* @param deltaS vector of rxn SS entropy changes
* units = J kmol-1 Kelvin-1
*/
virtual void getDeltaSSEntropy(doublereal* deltaS);
//@}
/**
* @name Species Production Rates
*/
//@{
//@}
/**
* @name Species Production Rates
*/
//@{
//! Returns the Species creation rates [kmol/m^2/s].
/*!
* Return the species
* creation rates in array cdot, which must be
* dimensioned at least as large as the total number of
* species in all phases of the kinetics
* model
*
* @param cdot Vector containing creation rates.
* length = m_kk. units = kmol/m^2/s
*/
virtual void getCreationRates(doublereal* cdot);
//! Returns the Species creation rates [kmol/m^2/s].
/*!
* Return the species
* creation rates in array cdot, which must be
* dimensioned at least as large as the total number of
* species in all phases of the kinetics
* model
*
* @param cdot Vector containing creation rates.
* length = m_kk. units = kmol/m^2/s
*/
virtual void getCreationRates(doublereal* cdot);
//! Return the Species destruction rates [kmol/m^2/s].
/*!
* Return the species destruction rates in array ddot, which must be
* dimensioned at least as large as the total number of
* species in all phases of the kinetics model
*
* @param ddot Vector containing destruction rates.
* length = m_kk. units = kmol/m^2/s
*/
virtual void getDestructionRates(doublereal* ddot);
//! Return the Species destruction rates [kmol/m^2/s].
/*!
* Return the species destruction rates in array ddot, which must be
* dimensioned at least as large as the total number of
* species in all phases of the kinetics model
*
* @param ddot Vector containing destruction rates.
* length = m_kk. units = kmol/m^2/s
*/
virtual void getDestructionRates(doublereal* ddot);
//! Return the species net production rates [kmol/m^2/s].
/*!
* Species net production rates [kmol/m^2/s]. Return the species
* net production rates (creation - destruction) in array
* wdot, which must be dimensioned at least as large as the
* total number of species in all phases of the kinetics
* model
*
* @param net Vector of species production rates.
* units kmol m-d s-1, where d is dimension.
*/
virtual void getNetProductionRates(doublereal* net);
//! Return the species net production rates [kmol/m^2/s].
/*!
* Species net production rates [kmol/m^2/s]. Return the species
* net production rates (creation - destruction) in array
* wdot, which must be dimensioned at least as large as the
* total number of species in all phases of the kinetics
* model
*
* @param net Vector of species production rates.
* units kmol m-d s-1, where d is dimension.
*/
virtual void getNetProductionRates(doublereal* net);
//@}
/**
* @name Reaction Mechanism Informational Query Routines
*/
//@{
//@}
/**
* @name Reaction Mechanism Informational Query Routines
*/
//@{
/**
* Stoichiometric coefficient of species k as a reactant in
* reaction i.
*/
virtual doublereal reactantStoichCoeff(int k, int i) const {
return m_rrxn[k][i];
}
/**
* Stoichiometric coefficient of species k as a reactant in
* reaction i.
*/
virtual doublereal reactantStoichCoeff(int k, int i) const {
return m_rrxn[k][i];
}
/**
* Stoichiometric coefficient of species k as a product in
* reaction i.
*/
virtual doublereal productStoichCoeff(int k, int i) const {
return m_prxn[k][i];
}
/**
* Stoichiometric coefficient of species k as a product in
* reaction i.
*/
virtual doublereal productStoichCoeff(int k, int i) const {
return m_prxn[k][i];
}
/**
* Flag specifying the type of reaction. The legal values and
* their meaning are specific to the particular kinetics
* manager.
*/
virtual int reactionType(int i) const {
return m_index[i].first;
}
/**
* Flag specifying the type of reaction. The legal values and
* their meaning are specific to the particular kinetics
* manager.
*/
virtual int reactionType(int i) const {
return m_index[i].first;
}
//! Return the charge transfer rxn Beta parameter for the ith reaction
/*!
* Returns the beta parameter for a charge transfer reaction. This
* parameter is not important for non-charge transfer reactions.
* Note, the parameter defaults to zero. However, a value of 0.5
* should be supplied for every charge transfer reaction if
* no information is known, as a value of 0.5 pertains to a
* symmetric transition state. The value can vary between 0 to 1.
*
*
* @param irxn Reaction number in the kinetics mechanism
*
* @return
* Beta parameter. This defaults to zero, even for charge transfer
* reactions.
*/
doublereal electrochem_beta(int irxn) const;
//! Return the charge transfer rxn Beta parameter for the ith reaction
/*!
* Returns the beta parameter for a charge transfer reaction. This
* parameter is not important for non-charge transfer reactions.
* Note, the parameter defaults to zero. However, a value of 0.5
* should be supplied for every charge transfer reaction if
* no information is known, as a value of 0.5 pertains to a
* symmetric transition state. The value can vary between 0 to 1.
*
*
* @param irxn Reaction number in the kinetics mechanism
*
* @return
* Beta parameter. This defaults to zero, even for charge transfer
* reactions.
*/
doublereal electrochem_beta(int irxn) const;
/**
* True if reaction i has been declared to be reversible. If
* isReversible(i) is false, then the reverse rate of progress
* for reaction i is always zero.
*/
virtual bool isReversible(int i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) return true;
else return false;
}
/**
* True if reaction i has been declared to be reversible. If
* isReversible(i) is false, then the reverse rate of progress
* for reaction i is always zero.
*/
virtual bool isReversible(int i) {
if (std::find(m_revindex.begin(), m_revindex.end(), i)
< m_revindex.end()) return true;
else return false;
}
/**
* Return a string representing the reaction.
*/
virtual std::string reactionString(int i) const {
return m_rxneqn[i];
}
/**
* Return a string representing the reaction.
*/
virtual std::string reactionString(int i) const {
return m_rxneqn[i];
}
virtual void getFwdRateConstants(doublereal* kfwd);
virtual void getRevRateConstants(doublereal* krev,
bool doIrreversible = false);
virtual void getFwdRateConstants(doublereal* kfwd);
virtual void getRevRateConstants(doublereal* krev,
bool doIrreversible = false);
virtual void getActivationEnergies(doublereal *E);
virtual void getActivationEnergies(doublereal *E);
//@}
/**
* @name Reaction Mechanism Construction
*/
//@{
//@}
/**
* @name Reaction Mechanism Construction
*/
//@{
//! Prepare the class for the addition of reactions.
/*!
* This function must be called after instantiation of the class, but before
* any reactions are actually added to the mechanism.
* This function calculates m_kk the number of species in all
* phases participating in the reaction mechanism. We don't know
* m_kk previously, before all phases have been added.
*/
virtual void init();
//! Prepare the class for the addition of reactions.
/*!
* This function must be called after instantiation of the class, but before
* any reactions are actually added to the mechanism.
* This function calculates m_kk the number of species in all
* phases participating in the reaction mechanism. We don't know
* m_kk previously, before all phases have been added.
*/
virtual void init();
//! Add a single reaction to the mechanism.
/*!
* @param r Reference to a ReactionData object containing all of
* the info needed to describe the reaction.
*/
virtual void addReaction(const ReactionData& r);
//! Add a single reaction to the mechanism.
/*!
* @param r Reference to a ReactionData object containing all of
* the info needed to describe the reaction.
*/
virtual void addReaction(const ReactionData& r);
//! Finish adding reactions and prepare for use.
/*!
* This function
* must be called after all reactions are entered into the mechanism
* and before the mechanism is used to calculate reaction rates.
*/
virtual void finalize();
//! Finish adding reactions and prepare for use.
/*!
* This function
* must be called after all reactions are entered into the mechanism
* and before the mechanism is used to calculate reaction rates.
*/
virtual void finalize();
virtual bool ready() const;
virtual bool ready() const;
void updateROP();
void updateROP();
void _update_rates_T();
void _update_rates_phi();
void _update_rates_C();
void _update_rates_T();
void _update_rates_phi();
void _update_rates_C();
//! Advance the surface coverages in time
/*!
* This method carries out a time-accurate advancement of the
* surface coverages for a specified amount of time.
*
* \f[
* \dot {\theta}_k = \dot s_k (\sigma_k / s_0)
* \f]
*
*
* @param tstep Time value to advance the surface coverages
*/
void advanceCoverages(doublereal tstep);
//! Advance the surface coverages in time
/*!
* This method carries out a time-accurate advancement of the
* surface coverages for a specified amount of time.
*
* \f[
* \dot {\theta}_k = \dot s_k (\sigma_k / s_0)
* \f]
*
*
* @param tstep Time value to advance the surface coverages
*/
void advanceCoverages(doublereal tstep);
//! Solve for the pseudo steady-state of the surface problem
/*!
* Solve for the steady state of the surface problem.
* This is the same thing as the advanceCoverages() function,
* but at infinite times.
*
* Note, a direct solve is carried out under the hood here,
* to reduce the computational time.
*/
void solvePseudoSteadyStateProblem();
//! Solve for the pseudo steady-state of the surface problem
/*!
* Solve for the steady state of the surface problem.
* This is the same thing as the advanceCoverages() function,
* but at infinite times.
*
* Note, a direct solve is carried out under the hood here,
* to reduce the computational time.
*
* @param ifuncOverride 4 values are possible
* 1 SFLUX_INITIALIZE
* 2 SFLUX_RESIDUAL
* 3 SFLUX_JACOBIAN
* 4 SFLUX_TRANSIENT
* The default is -1, which means that the program
* will decide.
* @param timeScaleOverride When a psuedo transient is
* selected this value can be used to override
* the default time scale for integration which
* is one.
* When SFLUX_TRANSIENT is used, this is equal to the
* time over which the equations are integrated.
* When SFLUX_INITIALIZE is used, this is equal to the
* time used in the initial transient algorithm,
* before the equation system is solved directly.
*/
void solvePseudoSteadyStateProblem(int ifuncOverride = -1,
doublereal timeScaleOverride = 1.0);
void checkPartialEquil();
void setIOFlag(int ioFlag);
void checkPartialEquil();
//! Temporary work vector of length m_kk
vector_fp m_grt;
//! Temporary work vector of length m_kk
vector_fp m_grt;
protected:
protected:
//! m_kk is the number of species in all of the phases
//! that participate in this kinetics mechanism.
int m_kk;
//! m_kk is the number of species in all of the phases
//! that participate in this kinetics mechanism.
int m_kk;
//! List of reactions numbers which are reversible reactions
/*!
* This is a vector of reaction numbers. Each reaction
* in the list is reversible.
* Length = number of reversible reactions
*/
vector_int m_revindex;
//! List of reactions numbers which are reversible reactions
/*!
* This is a vector of reaction numbers. Each reaction
* in the list is reversible.
* Length = number of reversible reactions
*/
vector_int m_revindex;
Rate1<SurfaceArrhenius> m_rates;
bool m_redo_rates;
Rate1<SurfaceArrhenius> m_rates;
bool m_redo_rates;
/**
* Vector of information about reactions in the
* mechanism.
* The key is the reaction index (0 < i < m_ii).
* The first pair is the reactionType of the reaction.
* The second pair is ...
*/
mutable std::map<int, std::pair<int, int> > m_index;
/**
* Vector of information about reactions in the
* mechanism.
* The key is the reaction index (0 < i < m_ii).
* The first pair is the reactionType of the reaction.
* The second pair is ...
*/
mutable std::map<int, std::pair<int, int> > m_index;
//! Vector of irreversible reaction numbers
/*!
* vector containing the reaction numbers of irreversible
* reactions.
*/
std::vector<int> m_irrev;
//! Vector of irreversible reaction numbers
/*!
* vector containing the reaction numbers of irreversible
* reactions.
*/
std::vector<int> m_irrev;
//! Stoichiometric manager for the reaction mechanism
/*!
* This is the manager for the kinetics mechanism that
* handles turning reaction extents into species
* production rates and also handles turning thermo
* properties into reaction thermo properties.
*/
ReactionStoichMgr m_rxnstoich;
//! Stoichiometric manager for the reaction mechanism
/*!
* This is the manager for the kinetics mechanism that
* handles turning reaction extents into species
* production rates and also handles turning thermo
* properties into reaction thermo properties.
*/
ReactionStoichMgr m_rxnstoich;
//! Number of irreversible reactions in the mechanism
int m_nirrev;
//! Number of irreversible reactions in the mechanism
int m_nirrev;
//! Number of reversible reactions in the mechanism
int m_nrev;
//! Number of reversible reactions in the mechanism
int m_nrev;
//! m_rrxn is a vector of maps, containing the reactant
//! stochiometric coefficient information
/*!
* m_rrxn has a length
* equal to the total number of species in the kinetics
* object. For each species, there exists a map, with the
* reaction number being the key, and the
* reactant stoichiometric coefficient for the species being the value.
* HKM -> mutable because search sometimes creates extra
* entries. To be fixed in future...
*/
mutable std::vector<std::map<int, doublereal> > m_rrxn;
//! m_rrxn is a vector of maps, containing the reactant
//! stochiometric coefficient information
/*!
* m_rrxn has a length
* equal to the total number of species in the kinetics
* object. For each species, there exists a map, with the
* reaction number being the key, and the
* reactant stoichiometric coefficient for the species being the value.
* HKM -> mutable because search sometimes creates extra
* entries. To be fixed in future...
*/
mutable std::vector<std::map<int, doublereal> > m_rrxn;
//! m_prxn is a vector of maps, containing the reactant
//! stochiometric coefficient information
/**
* m_prxn is a vector of maps. m_prxn has a length
* equal to the total number of species in the kinetics
* object. For each species, there exists a map, with the
* reaction number being the key, and the
* product stoichiometric coefficient for the species being the value.
*/
mutable std::vector<std::map<int, doublereal> > m_prxn;
//! String expression for each rxn
/*!
* Vector of strings of length m_ii, the number of
* reactions, containing the
* string expressions for each reaction
* (e.g., reactants <=> product1 + product2)
*/
std::vector<std::string> m_rxneqn;
//! m_prxn is a vector of maps, containing the reactant
//! stochiometric coefficient information
/**
* m_prxn is a vector of maps. m_prxn has a length
* equal to the total number of species in the kinetics
* object. For each species, there exists a map, with the
* reaction number being the key, and the
* product stoichiometric coefficient for the species being the value.
*/
mutable std::vector<std::map<int, doublereal> > m_prxn;
//! String expression for each rxn
/*!
* Vector of strings of length m_ii, the number of
* reactions, containing the
* string expressions for each reaction
* (e.g., reactants <=> product1 + product2)
*/
std::vector<std::string> m_rxneqn;
/**
* Temporary data storage used in calculating the rates of
* of reactions.
*/
InterfaceKineticsData* m_kdata;
/**
* Temporary data storage used in calculating the rates of
* of reactions.
*/
InterfaceKineticsData* m_kdata;
//! an array of generalized concentrations for each species
/*!
* An array of generalized concentrations
* \f$ C_k \f$ that are defined such that \f$ a_k = C_k /
* C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration/
* These generalized concentrations are used
* by this kinetics manager class to compute the forward and
* reverse rates of elementary reactions. The "units" for the
* concentrations of each phase depend upon the implementation
* of kinetics within that phase.
* The order of the species within the vector is based on
* the order of listed ThermoPhase objects in the class, and the
* order of the species within each ThermoPhase class.
*/
vector_fp m_conc;
//! an array of generalized concentrations for each species
/*!
* An array of generalized concentrations
* \f$ C_k \f$ that are defined such that \f$ a_k = C_k /
* C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration/
* These generalized concentrations are used
* by this kinetics manager class to compute the forward and
* reverse rates of elementary reactions. The "units" for the
* concentrations of each phase depend upon the implementation
* of kinetics within that phase.
* The order of the species within the vector is based on
* the order of listed ThermoPhase objects in the class, and the
* order of the species within each ThermoPhase class.
*/
vector_fp m_conc;
//! Vector of standard state chemical potentials
/*!
* This vector contains a temporary vector of
* standard state chemical potentials
* for all of the species in the kinetics object
*
* Length = m_k
* units = J/kmol
*/
vector_fp m_mu0;
//! Vector of standard state chemical potentials
/*!
* This vector contains a temporary vector of
* standard state chemical potentials
* for all of the species in the kinetics object
*
* Length = m_k
* units = J/kmol
*/
vector_fp m_mu0;
//! Vector of phase potentials
/*!
* Temporary vector containing the potential of each phase
* in the kinetics object
*
* length = number of phases
* units = Volts
*/
vector_fp m_phi;
//! Vector of phase potentials
/*!
* Temporary vector containing the potential of each phase
* in the kinetics object
*
* length = number of phases
* units = Volts
*/
vector_fp m_phi;
//! Vector of potential energies due to Voltages
/*!
* Length is the number of species in kinetics mech. It's
* used to store the potential energy due to the voltage.
*/
vector_fp m_pot;
//! Vector of potential energies due to Voltages
/*!
* Length is the number of species in kinetics mech. It's
* used to store the potential energy due to the voltage.
*/
vector_fp m_pot;
//! Vector temporary
/*!
* Length is number of reactions. it's used to store the
* voltage contribution to the activation energy.
*/
vector_fp m_rwork;
//! Vector temporary
/*!
* Length is number of reactions. it's used to store the
* voltage contribution to the activation energy.
*/
vector_fp m_rwork;
//! Vector of raw activation energies for the reactions
/*!
* units are in Kelvin
*/
vector_fp m_E;
//! Vector of raw activation energies for the reactions
/*!
* units are in Kelvin
*/
vector_fp m_E;
//! Pointer to the surface phase
SurfPhase* m_surf;
//! Pointer to the single surface phase
SurfPhase* m_surf;
//! Pointer to the surface solver
ImplicitSurfChem* m_integrator;
//! Pointer to the Implicit surface chemistry object
/*!
* Note this object is owned by this InterfaceKinetics
* object. It may only be used to solve this single
* InterfaceKinetics objects's surface problem uncoupled
* from other surface phases.
*/
ImplicitSurfChem* m_integrator;
vector_fp m_beta;
vector_int m_ctrxn;
vector_fp m_beta;
vector_int m_ctrxn;
int reactionNumber(){ return m_ii;}
int reactionNumber(){ return m_ii;}
void addElementaryReaction(const ReactionData& r);
void addGlobalReaction(const ReactionData& r);
void installReagents(const ReactionData& r);
void addElementaryReaction(const ReactionData& r);
void addGlobalReaction(const ReactionData& r);
void installReagents(const ReactionData& r);
void updateKc();
void updateKc();
//! Write values into m_index
/*!
* @param rxnNumber reaction number
* @param type reaction type
* @param loc location ??
*/
void registerReaction(int rxnNumber, int type, int loc) {
m_index[rxnNumber] = std::pair<int, int>(type, loc);
}
//! Write values into m_index
/*!
* @param rxnNumber reaction number
* @param type reaction type
* @param loc location ??
*/
void registerReaction(int rxnNumber, int type, int loc) {
m_index[rxnNumber] = std::pair<int, int>(type, loc);
}
void applyButlerVolmerCorrection(doublereal* kf);
void applyButlerVolmerCorrection(doublereal* kf);
//! boolean indicating whether mechanism has been finalized
bool m_finalized;
bool m_has_coverage_dependence;
bool m_has_electrochem_rxns;
//! boolean indicating whether mechanism has been finalized
bool m_finalized;
bool m_has_coverage_dependence;
bool m_has_electrochem_rxns;
private:
int m_ioFlag;
private:
};
};
}
#endif

6
Cantera/src/kinetics/Makefile.in
View file

@ -27,7 +27,7 @@ else
endif
PIC_FLAG=@PIC@
LOCAL_DEFS = -DDEBUG_SOLVESP
CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG)
# homogeneous kinetics
@ -35,12 +35,12 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG)
ifeq ($(do_kinetics),1)
KINETICS_OBJ=importKinetics.o GRI_30_Kinetics.o KineticsFactory.o \
GasKinetics.o \
FalloffFactory.o ReactionStoichMgr.o Kinetics.o
FalloffFactory.o ReactionStoichMgr.o Kinetics.o solveSP.o
KINETICS_H = importKinetics.h GRI_30_Kinetics.h KineticsFactory.h \
Kinetics.h GasKinetics.h \
FalloffFactory.h ReactionStoichMgr.h reaction_defs.h \
FalloffMgr.h ThirdBodyMgr.h RateCoeffMgr.h ReactionData.h \
RxnRates.h Enhanced3BConc.h StoichManager.h
RxnRates.h Enhanced3BConc.h StoichManager.h solveSP.h
KINETICS = $(KINETICS_OBJ) $(KINETICS_H)
endif

1445
Cantera/src/kinetics/solveSP.cpp Normal file
View file

File diff suppressed because it is too large Load diff

686
Cantera/src/kinetics/solveSP.h Normal file
View file

@ -0,0 +1,686 @@
/**
* @file solveSP.h
* Header file for implicit surface problem solver
* (see \ref kinetics and class \link Cantera::solveSP solveSP\endlink).
*/
/*
* $Id$
*/
/*
* Copywrite 2004 Sandia Corporation. Under the terms of Contract
* DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
* retains certain rights in this software.
* See file License.txt for licensing information.
*/
#ifndef SOLVESP_H
#define SOLVESP_H
#include "ImplicitSurfChem.h"
#include "InterfaceKinetics.h"
#include <vector>
#include "Array.h"
//! Solution Methods
/*!
* Flag to specify the solution method
*
* 1: SFLUX_INITIALIZE = This assumes that the initial guess supplied to the
* routine is far from the correct one. Substantial
* work plus transient time-stepping is to be expected
* to find a solution.
* 2: SFLUX_RESIDUAL = Need to solve the surface problem in order to
* calculate the surface fluxes of gas-phase species.
* (Can expect a moderate change in the solution
* vector -> try to solve the system by direct
* methods
* with no damping first -> then, try time-stepping
* if the first method fails)
* A "time_scale" supplied here is used in the
* algorithm to determine when to shut off
* time-stepping.
* 3: SFLUX_JACOBIAN = Calculation of the surface problem is due to the
* need for a numerical jacobian for the gas-problem.
* The solution is expected to be very close to the
* initial guess, and accuracy is needed.
* 4: SFLUX_TRANSIENT = The transient calculation is performed here for an
* amount of time specified by "time_scale". It is
* not garraunted to be time-accurate - just stable
* and fairly fast. The solution after del_t time is
* returned, whether it's converged to a steady
* state or not.
*/
const int SFLUX_INITIALIZE = 1;
const int SFLUX_RESIDUAL = 2;
const int SFLUX_JACOBIAN = 3;
const int SFLUX_TRANSIENT = 4;
/*
* bulkFunc: Functionality expected from the bulk phase. This changes the
* equations that will be used to solve for the bulk mole
* fractions.
* 1: BULK_DEPOSITION = deposition of a bulk phase is to be expected.
* Bulk mole fractions are determined from ratios of
* growth rates of bulk species.
* 2: BULK_ETCH = Etching of a bulk phase is to be expected.
* Bulk mole fractions are assumed constant, and given
* by the initial conditions. This is also used
whenever the condensed phase is part of the larger
solution.
*/
const int BULK_DEPOSITION = 1;
const int BULK_ETCH = 2;
namespace Cantera {
class InterfaceKinetics;
//! Method to solve a pseudo steady state surface problem
/*!
* The following class handles solving the surface problem.
* The calculation uses Newton's method to
* obtain the surface fractions of the surface and bulk species by
* requiring that the
* surface species production rate = 0 and that the either the
* bulk fractions are proportional to their production rates
* or they are constants.
*
* Currently, the bulk mole fractions are treated as constants.
* Implementation of their being added to the unknown solution
* vector is delayed.
*
* Lets introduce the unknown vector for the "surface
* problem". The surface problem is defined as the evaluation of the surface
* site fractions for multiple surface phases.
* The unknown vector will consist of the vector of surface concentrations for each
* species in each surface vector. Species are grouped first by their surface phases
*
* C_i_j = Concentration of the ith species in the jth surface phase
* Nj = number of surface species in the jth surface phase
*
* The unknown solution vector is defined as follows:
*
* kindexSP
* ----------------------------
* C_0_0 0
* C_1_0 1
* C_2_0 2
* . . . ...
* C_N0-1_0 N0-1
* C_0_1 N0
* C_1_1 N0+1
* C_2_1 N0+2
* . . . ...
* C_N1-1_1 NO+N1-1
*
*
* Note there are a couple of different types of species indecices
* floating around in the formulation of this object.
*
* kindexSP This is the species index in the contiguous vector of unknowns
* for the surface problem.
*
* Note, in the future, BULK_DEPOSITION systems will be added, and the solveSP unknown
* vector will get more complicated. It will include the mole fraction and growth rates
* of specified bulk phases
*
* Indecises which relate to individual kinetics objects use the suffix KSI (kinetics
* species index).
*
*
* Solution Method
*
* This routine is typically used within a residual calculation in a large code.
* It's typically invoked millions of times for large calculations, and it must
* work every time. Therefore, requirements demand that it be robust but also
* efficient.
*
* The solution methodology is largely determined by the <TT>ifunc<\TT> parameter,
* that is input to the solution object. This parameter may have the following
* 4 values:
*
*
* 1: SFLUX_INITIALIZE = This assumes that the initial guess supplied to the
* routine is far from the correct one. Substantial
* work plus transient time-stepping is to be expected
* to find a solution.
*
* 2: SFLUX_RESIDUAL = Need to solve the surface problem in order to
* calculate the surface fluxes of gas-phase species.
* (Can expect a moderate change in the solution
* vector -> try to solve the system by direct methods
* with no damping first -> then, try time-stepping
* if the first method fails)
* A "time_scale" supplied here is used in the
* algorithm to determine when to shut off
* time-stepping.
*
* 3: SFLUX_JACOBIAN = Calculation of the surface problem is due to the
* need for a numerical jacobian for the gas-problem.
* The solution is expected to be very close to the
* initial guess, and extra accuracy is needed because
* solution variables have been delta'd from
* nominal values to create jacobian entries.
*
* 4: SFLUX_TRANSIENT = The transient calculation is performed here for an
* amount of time specified by "time_scale". It is
* not garraunted to be time-accurate - just stable
* and fairly fast. The solution after del_t time is
* returned, whether it's converged to a steady
* state or not. This is a poor man's time stepping
* algorithm.
*
* Psuedo time stepping algorithm:
* The time step is determined from sdot[], so so that the time step
* doesn't ever change the value of a variable by more than 100%.
*
* This algorithm does use a damped Newton's method to relax the equations.
* Damping is based on a "delta damping" technique. The solution unknowns
* are not allowed to vary too much between iterations.
*
*
* EXTRA_ACCURACY:A constant that is the ratio of the required update norm in
* this Newton iteration compared to that in the nonlinear solver.
* A value of 0.1 is used so surface species are safely overconverged.
*
* Functions called:
*----------------------------------------------------------------------------
*
* ct_dgetrf -- First half of LAPACK direct solve of a full Matrix
*
* ct_dgetrs -- Second half of LAPACK direct solve of a full matrix. Returns
* solution vector in the right-hand-side vector, resid.
*
*----------------------------------------------------------------------------
*
*/
class solveSP {
public:
//! Constructor for the object
/*!
* @param surfChemPtr Pointer to the ImplicitSurfChem object that
* defines the surface problem to be solved.
*
* @param bulkFunc Integer representing how the bulk phases
* should be handled. Currently, only the
* default value of BULK_ETCH is supported.
*/
solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc = BULK_ETCH);
//! Destructor. Deletes the integrator.
~solveSP();
private:
//! Unimplemented private copy constructor
solveSP(const solveSP &right);
//! Unimplemented private assignment operator
solveSP& operator=(const solveSP &right);
public:
//! Main routine that actually calculates the pseudo steady state
//! of the surface problem
/*!
* The actual converged solution is returned as part of the
* internal state of the InterfaceKinetics objects.
*
* @param ifunc Determines the type of solution algorithm to be
* used. Possible values are SFLUX_INITIALIZE ,
* SFLUX_RESIDUAL SFLUX_JACOBIAN SFLUX_TRANSIENT .
*
* @param time_scale Time over which to integrate the surface equations,
* where applicable
*
* @param TKelvin Temperature (kelvin)
*
* @param PGas Pressure (pascals)
*
* @param reltol Relative tolerance to use
* @param abstol absolute tolerance.
*
* @return Returns 1 if the surface problem is successfully solved.
* Returns -1 if the surface problem wasn't solved successfully.
* Note the actual converged solution is returned as part of the
* internal state of the InterfaceKinetics objects.
*/
int solveSurfProb(int ifunc, double time_scale, double TKelvin,
double PGas, double reltol, double abstol);
private:
//! Printing routine that gets called at the start of every
//! invocation
void print_header(int ioflag, int ifunc, double time_scale,
int damping, double reltol, double abstol,
double TKelvin, double PGas, double netProdRate[],
double XMolKinSpecies[]);
#ifdef DEBUG_SOLVESP
void printResJac(int ioflag, int neq, const Array2D &Jac,
double resid[], double wtResid[], double norm);
#endif
//! Printing routine that gets called after every iteration
void printIteration(int ioflag, double damp, int label_d, int label_t,
double inv_t, double t_real, int iter,
double update_norm, double resid_norm,
double netProdRate[], double CSolnSP[],
double resid[], double XMolSolnSP[],
double wtSpecies[], int dim, bool do_time);
//! Print a summary of the solution
/*!
*
*/
void printFinal(int ioflag, double damp, int label_d, int label_t,
double inv_t, double t_real, int iter,
double update_norm, double resid_norm,
double netProdRateKinSpecies[], const double CSolnSP[],
const double resid[], double XMolSolnSP[],
const double wtSpecies[], const double wtRes[],
int dim, bool do_time,
double TKelvin, double PGas);
//! Calculate a conservative delta T to use in a pseudo-steady state
//! algorithm
/*!
* This routine calculates a pretty conservative 1/del_t based
* on MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees
* diagonal dominance.
*
* Small surface fractions are allowed to intervene in the del_t
* determination, no matter how small. This may be changed.
* Now minimum changed to 1.0e-12,
*
* Maximum time step set to time_scale.
*
* @param netProdRateSolnSP Output variable. Net production rate
* of all of the species in the solution vector.
* @param XMolSolnSP output variable.
* Mole fraction of all of the species in the solution vector
* @param label Output variable. Pointer to the value of the
* species index (kindexSP) that is controlling
* the time step
* @param label_old Output variable. Pointer to the value of the
* species index (kindexSP) that controlled
* the time step at the previous iteration
* @param label_factor Output variable. Pointer to the current
* factor that is used to indicate the same species
* is controlling the time step.
*
* @param ioflag Level of the output requested.
*
* @return Returns the 1. / delta T to be used on the next step
*/
double calc_t(double netProdRateSolnSP[], double XMolSolnSP[],
int *label, int *label_old,
double *label_factor, int ioflag);
//! Calculate the solution and residual weights
/*!
* @param wtSpecies Weights to use for the soln unknowns. These
* are in concentration units
* @param wtResid Weights to sue for the residual unknowns.
*
* @param Jac Jacobian. Row sum scaling is used for the Jacobian
* @param CSolnSP Solution vector for the surface problem
* @param abstol Absolute error tolerance
* @param reltol Relative error tolerance
*/
void calcWeights(double wtSpecies[], double wtResid[],
const Array2D &Jac, const double CSolnSP[],
const double abstol, const double reltol);
#ifdef DEBUG_SOLVESP
//! Utility routine to print a header for high lvls of debugging
/*!
* @param ioflag Lvl of debugging
* @param damp lvl of damping
* @param inv_t Inverse of the value of delta T
* @param t_real Value of the time
* @param iter Interation number
* @param do_time boolean indicating whether time stepping is taking
* place
*/
void printIterationHeader(int ioflag, double damp,
double inv_t, double t_real, int iter,
bool do_time);
#endif
/**
* Update the surface states of the surface phases.
*/
void updateState(const double* cSurfSpec);
//! Update mole fraction vector consisting of unknowns in surface problem
/*!
* @param XMolSolnSP Vector of mole fractions for the unknowns in the
* surface problem.
*/
void updateMFSolnSP(double *XMolSolnSP);
//! Update the mole fraction vector for a specific kinetic species vector
//! corresponding to one InterfaceKinetics object
/*!
* @param XMolKinSp Mole fraction vector corresponding to a particular
* kinetic species for a single InterfaceKinetics Object
* This is a vector over all the species in all of the
* phases in the InterfaceKinetics object
* @param isp ID of the InterfaceKinetics Object.
*/
void updateMFKinSpecies(double *XMolKinSp, int isp);
//! Update the vector that keeps track of the largest species in each
//! surface phase.
/*!
* @param CsolnSP Vector of the current values of the surface concentrations
* in all of the surface species.
*/
void evalSurfLarge(const double *CSolnSP);
//! Main Function evalulation
/*!
*
* @param resid output Vector of residuals, length = m_neq
* @param CSolnSP Vector of species concentrations, unknowns in the
* problem, length = m_neq
* @param CSolnSPOld Old Vector of species concentrations, unknowns in the
* problem, length = m_neq
* @param do_time Calculate a time dependent residual
* @param deltaT Delta time for time dependent problem.
*/
void fun_eval(double* resid, const double *CSolnSP,
const double *CSolnOldSP, const bool do_time, const double deltaT);
//! Main routine that calculates the current residual and Jacobian
/*!
* @param JacCol Vector of pointers to the tops of columns of the
* Jacobian to be evalulated.
* @param resid output Vector of residuals, length = m_neq
* @param CSolnSP Vector of species concentrations, unknowns in the
* problem, length = m_neq. These are tweaked in order
* to derive the columns of the jacobian.
* @param CSolnSPOld Old Vector of species concentrations, unknowns in the
* problem, length = m_neq
* @param do_time Calculate a time dependent residual
* @param deltaT Delta time for time dependent problem.
*/
void resjac_eval(std::vector<double*>& JacCol, double* resid,
double *CSolnSP,
const double *CSolnSPOld, const bool do_time,
const double deltaT);
//! Pointer to the manager of the implicit surface chemistry
//! problem
/*!
* This object actually calls the current object. Thus, we are
* providing a loop-back functionality here.
*/
ImplicitSurfChem *m_SurfChemPtr;
//! Vector of interface kinetics objects
/*!
* Each of these is associated with one and only one surface phase.
*/
std::vector<InterfaceKinetics*> &m_objects;
//! Total number of equations to solve in the implicit problem.
/*!
* Note, this can be zero, and frequently is
*/
int m_neq;
//! This variable determines how the bulk phases are to be handled
/*!
* = BULK_ETCH (default) The concentrations of the bulk phases are
* considered constant, just as the gas phase is.
* They are not part of the solution vector.
* = BULK_DEPOSITION =
* We solve here for the composition of the bulk
* phases by calculating a growth rate. The equations
* for the species in the bulk phases are
* unknowns in this calculation.
*/
int m_bulkFunc;
//! Number of surface phases in the surface problem
/*!
* This number is equal to the number of InterfaceKinetics objects
* in the problem. (until further noted)
*/
int m_numSurfPhases;
//! Total number of surface species in all surface phases.
/*!
* This is also the number of equations to solve for m_mode=0 system
* It's equal to the sum of the number of species in each of the
* m_numSurfPhases.
*/
int m_numTotSurfSpecies;
//! Mapping between the surface phases and the InterfaceKinetics objects
/*!
* Currently this is defined to be a 1-1 mapping (and probably assumed
* in some places)
* m_surfKinObjID[i] = i
*/
std::vector<int> m_indexKinObjSurfPhase;
//! Vector of length number of surface phases containing
//! the number of surface species in each phase
/*!
* Length is equal to the number of surface phases, m_numSurfPhases
*/
std::vector<int> m_nSpeciesSurfPhase;
//! Vector of surface phase pointers
/*!
* This is created during the constructor
* Length is equal to the number of surface phases, m_numSurfPhases
*/
std::vector<SurfPhase *> m_ptrsSurfPhase;
//! Index of the start of the unknowns for each solution phase
/*!
* i_eqn = m_eqnIndexStartPhase[isp]
*
* isp is the phase id in the list of phases solved by the
* surface problem.
*
* i_eqn is the equation number of the first unknown in the
* solution vector corresponding to isp'th phase.
*/
std::vector<int> m_eqnIndexStartSolnPhase;
//! Phase ID in the InterfaceKinetics object of the surface phase
/*!
* For each surface phase, this lists the PhaseId of the
* surface phase in the corresponding InterfaceKinetics object
*
* Length is equal to m_numSurfPhases
*/
std::vector<int> m_kinObjPhaseIDSurfPhase;
//! Total number of volumetric condensed phases included in the steady state
//! problem handled by this routine.
/*!
* This is equal to or less
* than the total number of volumetric phases in all of the InterfaceKinetics
* objects. We usually do not include bulk phases. Bulk phases
* are only included in the calculation when their domain isn't included
* in the underlying continuum model conservation equation system.
*
* This is equal to 0, for the time being
*/
int m_numBulkPhasesSS;
//! Vector of number of species in the m_numBulkPhases phases.
/*!
* Length is number of bulk phases
*/
std::vector<int> m_numBulkSpecies;
//std::vector<int> m_bulkKinObjID;
//std::vector<int> m_bulkKinObjPhaseID;
//! Total number of species in all bulk phases.
/*!
* This is also the number of bulk equations to solve when bulk
* equation solving is turned on.
*/
int m_numTotBulkSpeciesSS;
//! Vector of bulk phase pointers, length is equal to m_numBulkPhases.
/*!
*
*/
std::vector<ThermoPhase *> m_bulkPhasePtrs;
//! Index between the equation index and the position in the
//! kinetic species array for the appropriate kinetics
//! operator
/*!
* Length = m_neq.
*
* ksp = m_kinSpecIndex[ieq]
* ksp is the kinetic species index for the ieq'th equation.
*/
std::vector<int> m_kinSpecIndex;
//! Index between the equation index and the index of the
//! InterfaceKinetics object
/*!
* Length m_neq
*/
std::vector<int> m_kinObjIndex;
//! Vector containing the indecies of the largest species
//! in each surface phase
/*!
* k = m_spSurfLarge[i]
* where
* k is the local species index, i.e.,
* it varies from 0 num species in phase-1
* i is the surface phase index in the problem
*
* length is equal to m_numSurfPhases
*/
std::vector<int> m_spSurfLarge;
//! m_atol is the absolute tolerance in real units.
/*!
* units are (kmol/m2)
*/
double m_atol;
//! m_rtol is the relative error tolerance.
double m_rtol;
//! maximum value of the time step
/*!
* units = seconds
*/
doublereal m_maxstep;
//! Maximum number of species in any single kinetics operator
//! -> also maxed wrt the total # of solution species
int m_maxTotSpecies;
//! Temporary vector with length equal to max m_maxTotSpecies
vector_fp m_netProductionRatesSave;
//! Temporary vector with length equal to max m_maxTotSpecies
vector_fp m_numEqn1;
//! Temporary vector with length equal to max m_maxTotSpecies
vector_fp m_numEqn2;
//! Temporary vector with length equal to max m_maxTotSpecies
vector_fp m_CSolnSave;
//! Solution vector
/*!
* length MAX(1, m_neq)
*/
vector_fp m_CSolnSP;
//! Saved inital solution vector
/*!
* length MAX(1, m_neq)
*/
vector_fp m_CSolnSPInit;
//! Saved solution vector at the old time step
/*!
* length MAX(1, m_neq)
*/
vector_fp m_CSolnSPOld;
//! Weights for the residual norm calculation
/*!
* length MAX(1, m_neq)
*/
vector_fp m_wtResid;
//! Weights for the species concentrations norm calculation
/*!
* length MAX(1, m_neq)
*/
vector_fp m_wtSpecies;
//! Residual for the surface problem
/*!
* The residual vector of length "dim" that, that has the value
* of "sdot" for surface species. The residuals for the bulk
* species are a function of the sdots for all species in the bulk
* phase. The last residual of each phase enforces {Sum(fractions)
* = 1}. After linear solve (dgetrf_ & dgetrs_), resid holds the
* update vector.
*
* length MAX(1, m_neq)
*/
vector_fp m_resid;
//! Vector of mole fractions
/*!
*length m_maxTotSpecies
*/
vector_fp m_XMolKinSpecies;
//! pivots
/*!
* length MAX(1, m_neq)
*/
vector_int m_ipiv;
//! Vector of pointers to the top of the columns of the
//! jacobians
/*!
* The "dim" by "dim" computed Jacobian matrix for the
* local Newton's method.
*/
std::vector<double *> m_JacCol;
//! Jacobian
/*!
* m_neq by m_neq computed Jacobian matrix for the
* local Newton's method.
*/
Array2D m_Jac;
public:
int ioflag;
};
}
#endif