Clean up comments in Kinetics
This commit is contained in:
Ray Speth
parent
5bf3b0fc67
commit
75b0c71044
22 changed files with 643 additions and 876 deletions
|
|
@ -17,9 +17,9 @@ namespace Cantera
|
|||
{
|
||||
|
||||
/**
|
||||
* Kinetics manager for elementary gas-phase chemistry. This
|
||||
* kinetics manager implements standard mass-action reaction rate
|
||||
* expressions for low-density gases.
|
||||
* Kinetics manager for elementary gas-phase chemistry. This kinetics manager
|
||||
* implements standard mass-action reaction rate expressions for low-density
|
||||
* gases.
|
||||
* @ingroup kinetics
|
||||
*/
|
||||
class GasKinetics : public BulkKinetics
|
||||
|
|
|
|||
|
|
@ -56,10 +56,9 @@ class ImplicitSurfChem : public FuncEval
|
|||
public:
|
||||
//! Constructor for multiple surfaces.
|
||||
/*!
|
||||
* @param k Vector of pointers to InterfaceKinetics objects
|
||||
* Each object consists of a surface or an edge containing
|
||||
* internal degrees of freedom representing the concentration
|
||||
* of surface adsorbates.
|
||||
* @param k Vector of pointers to InterfaceKinetics objects Each object
|
||||
* consists of a surface or an edge containing internal degrees of
|
||||
* freedom representing the concentration of surface adsorbates.
|
||||
*/
|
||||
ImplicitSurfChem(std::vector<InterfaceKinetics*> k);
|
||||
|
||||
|
|
@ -101,15 +100,13 @@ public:
|
|||
*
|
||||
* @param ifuncOverride One of the values defined in @ref solvesp_methods.
|
||||
* The default is -1, which means that the program will decide.
|
||||
* @param timeScaleOverride When a pseudo 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.
|
||||
* @param timeScaleOverride When a pseudo 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);
|
||||
|
|
@ -138,9 +135,8 @@ public:
|
|||
*
|
||||
* @param t0 Initial time
|
||||
* @param leny Length of the solution vector
|
||||
* @param y Value of the solution vector to be used.
|
||||
* On output, this contains the initial value
|
||||
* of the solution.
|
||||
* @param y Value of the solution vector to be used. On output, this
|
||||
* contains the initial value of the solution.
|
||||
* @deprecated Use getState() instead. To be removed after Cantera 2.3.
|
||||
*/
|
||||
virtual void getInitialConditions(doublereal t0,
|
||||
|
|
@ -161,9 +157,9 @@ public:
|
|||
* 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.
|
||||
* @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;
|
||||
|
||||
|
|
@ -173,10 +169,9 @@ public:
|
|||
* 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.
|
||||
* @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);
|
||||
|
||||
|
|
@ -272,31 +267,27 @@ protected:
|
|||
*/
|
||||
int m_mediumSpeciesStart;
|
||||
/**
|
||||
* Index into the species vector of the kinetics manager,
|
||||
* pointing to the first species from the condensed phase
|
||||
* of the particles.
|
||||
* 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
|
||||
* 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.
|
||||
* Pointer to the helper method, Placid, which solves the surface problem.
|
||||
*/
|
||||
std::unique_ptr<solveSP> m_surfSolver;
|
||||
|
||||
//! If true, a common temperature and pressure for all
|
||||
//! surface and bulk phases associated with the surface problem
|
||||
//! is imposed
|
||||
//! 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.
|
||||
//! 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:
|
||||
|
|
|
|||
|
|
@ -20,39 +20,37 @@ namespace Cantera
|
|||
class SurfPhase;
|
||||
class ImplicitSurfChem;
|
||||
|
||||
//! A kinetics manager for heterogeneous reaction mechanisms. The
|
||||
//! reactions are assumed to occur at a 2D interface between two 3D phases.
|
||||
//! A kinetics manager for heterogeneous reaction mechanisms. The reactions are
|
||||
//! assumed to occur at a 2D interface between two 3D phases.
|
||||
/*!
|
||||
* There are some important additions to the behavior of the kinetics class
|
||||
* due to the presence of multiple phases and a heterogeneous interface. If
|
||||
* a reactant phase doesn't exists, i.e., has a mole number of zero, a
|
||||
* heterogeneous reaction can not proceed from reactants to products. Note it
|
||||
* could perhaps proceed from products to reactants if all of the product
|
||||
* phases exist.
|
||||
* There are some important additions to the behavior of the kinetics class due
|
||||
* to the presence of multiple phases and a heterogeneous interface. If a
|
||||
* reactant phase doesn't exists, i.e., has a mole number of zero, a
|
||||
* heterogeneous reaction can not proceed from reactants to products. Note it
|
||||
* could perhaps proceed from products to reactants if all of the product phases
|
||||
* exist.
|
||||
*
|
||||
* In order to make the determination of whether a phase exists or not
|
||||
* actually involves the specification of additional information to the
|
||||
* kinetics object., which heretofore has only had access to intrinsic field
|
||||
* information about the phases (i.e., temperature pressure, and mole
|
||||
* fraction).
|
||||
* In order to make the determination of whether a phase exists or not actually
|
||||
* involves the specification of additional information to the kinetics object.,
|
||||
* which heretofore has only had access to intrinsic field information about the
|
||||
* phases (i.e., temperature pressure, and mole fraction).
|
||||
*
|
||||
* The extrinsic specification of whether a phase exists or not must be
|
||||
* specified on top of the intrinsic calculation of the reaction rate. This
|
||||
* class carries a set of booleans indicating whether a phase in the
|
||||
* heterogeneous mechanism exists or not.
|
||||
* The extrinsic specification of whether a phase exists or not must be
|
||||
* specified on top of the intrinsic calculation of the reaction rate. This
|
||||
* class carries a set of booleans indicating whether a phase in the
|
||||
* heterogeneous mechanism exists or not.
|
||||
*
|
||||
* Additionally, the class carries a set of booleans around indicating
|
||||
* whether a product phase is stable or not. If a phase is not
|
||||
* thermodynamically stable, it may be the case that a particular reaction in
|
||||
* a heterogeneous mechanism will create a product species in the unstable
|
||||
* phase. However, other reactions in the mechanism will destruct that
|
||||
* species. This may cause oscillations in the formation of the unstable
|
||||
* phase from time step to time step within a ODE solver, in practice. In
|
||||
* order to avoid this situation, a set of booleans is tracked which sets the
|
||||
* stability of a phase. If a phase is deemed to be unstable, then species in
|
||||
* that phase will not be allowed to be birthed by the kinetics operator.
|
||||
* Nonexistent phases are deemed to be unstable by default, but this can be
|
||||
* changed.
|
||||
* Additionally, the class carries a set of booleans around indicating whether a
|
||||
* product phase is stable or not. If a phase is not thermodynamically stable,
|
||||
* it may be the case that a particular reaction in a heterogeneous mechanism
|
||||
* will create a product species in the unstable phase. However, other reactions
|
||||
* in the mechanism will destruct that species. This may cause oscillations in
|
||||
* the formation of the unstable phase from time step to time step within a ODE
|
||||
* solver, in practice. In order to avoid this situation, a set of booleans is
|
||||
* tracked which sets the stability of a phase. If a phase is deemed to be
|
||||
* unstable, then species in that phase will not be allowed to be birthed by the
|
||||
* kinetics operator. Nonexistent phases are deemed to be unstable by default,
|
||||
* but this can be changed.
|
||||
*
|
||||
* @ingroup chemkinetics
|
||||
*/
|
||||
|
|
@ -61,12 +59,11 @@ 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.
|
||||
* @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);
|
||||
|
||||
|
|
@ -102,20 +99,10 @@ public:
|
|||
*/
|
||||
virtual void getEquilibriumConstants(doublereal* kc);
|
||||
|
||||
/** values needed to convert from exchange current density to surface reaction rate.
|
||||
*/
|
||||
//! values needed to convert from exchange current density to surface
|
||||
//! reaction rate.
|
||||
void updateExchangeCurrentQuantities();
|
||||
|
||||
//! Return the vector of values for the reaction Gibbs free energy change.
|
||||
/*!
|
||||
* (virtual from Kinetics.h)
|
||||
* These values depend upon the concentration of the solution.
|
||||
*
|
||||
* units = J kmol-1
|
||||
*
|
||||
* @param deltaG Output vector of deltaG's for reactions Length:
|
||||
* nReactions(). If 0, this updates the internally stored values only.
|
||||
*/
|
||||
virtual void getDeltaGibbs(doublereal* deltaG);
|
||||
|
||||
virtual void getDeltaElectrochemPotentials(doublereal* deltaM);
|
||||
|
|
@ -245,10 +232,9 @@ public:
|
|||
//! Update properties that depend on the species mole fractions and/or
|
||||
//! concentration,
|
||||
/*!
|
||||
* This method fills out the array of generalized concentrations by
|
||||
* calling method getActivityConcentrations for each phase, which classes
|
||||
* representing phases should overload to return the appropriate
|
||||
* quantities.
|
||||
* This method fills out the array of generalized concentrations by calling
|
||||
* method getActivityConcentrations for each phase, which classes
|
||||
* representing phases should overload to return the appropriate quantities.
|
||||
*/
|
||||
void _update_rates_C();
|
||||
|
||||
|
|
@ -275,15 +261,13 @@ public:
|
|||
*
|
||||
* @param ifuncOverride One of the values defined in @ref solvesp_methods.
|
||||
* The default is -1, which means that the program will decide.
|
||||
* @param timeScaleOverride When a pseudo 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.
|
||||
* @param timeScaleOverride When a pseudo 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);
|
||||
|
|
@ -292,22 +276,25 @@ public:
|
|||
|
||||
void checkPartialEquil();
|
||||
|
||||
//! Update the standard state chemical potentials and species equilibrium constant entries
|
||||
//! Update the standard state chemical potentials and species equilibrium
|
||||
//! constant entries
|
||||
/*!
|
||||
* Virtual because it is overwritten when dealing with experimental open circuit voltage overrides
|
||||
* Virtual because it is overwritten when dealing with experimental open
|
||||
* circuit voltage overrides
|
||||
*/
|
||||
virtual void updateMu0();
|
||||
|
||||
//! Update the equilibrium constants and stored electrochemical potentials
|
||||
//! in molar units for all reversible reactions and for all species.
|
||||
/*!
|
||||
* Irreversible reactions have their equilibrium constant set
|
||||
* to zero. For reactions involving charged species the equilibrium
|
||||
* constant is adjusted according to the electrostatic potential.
|
||||
* Irreversible reactions have their equilibrium constant set
|
||||
* to zero. For reactions involving charged species the equilibrium
|
||||
* constant is adjusted according to the electrostatic potential.
|
||||
*/
|
||||
void updateKc();
|
||||
|
||||
//! Apply modifications for the forward reaction rate for interfacial charge transfer reactions
|
||||
//! Apply modifications for the forward reaction rate for interfacial charge
|
||||
//! transfer reactions
|
||||
/*!
|
||||
* For reactions that transfer charge across a potential difference,
|
||||
* the activation energies are modified by the potential difference.
|
||||
|
|
@ -319,16 +306,17 @@ public:
|
|||
void applyVoltageKfwdCorrection(doublereal* const kfwd);
|
||||
|
||||
//! When an electrode reaction rate is optionally specified in terms of its
|
||||
//! exchange current density, adjust kfwd to the standard reaction rate constant form and units.
|
||||
//! When the BV reaction types are used, keep the exchange current density form.
|
||||
//! exchange current density, adjust kfwd to the standard reaction rate
|
||||
//! constant form and units. When the BV reaction types are used, keep the
|
||||
//! exchange current density form.
|
||||
/*!
|
||||
* For a reaction rate constant that was given in units of Amps/m2 (exchange current
|
||||
* density formulation with iECDFormulation == true), convert the rate to
|
||||
* kmoles/m2/s.
|
||||
* For a reaction rate constant that was given in units of Amps/m2
|
||||
* (exchange current density formulation with iECDFormulation == true),
|
||||
* convert the rate to kmoles/m2/s.
|
||||
*
|
||||
* For a reaction rate constant that was given in units of kmol/m2/sec when the
|
||||
* reaction type is a Butler-Volmer form, convert it to exchange current density
|
||||
* form (amps/m2).
|
||||
* For a reaction rate constant that was given in units of kmol/m2/sec when
|
||||
* the reaction type is a Butler-Volmer form, convert it to exchange
|
||||
* current density form (amps/m2).
|
||||
*
|
||||
* @param kfwd Vector of forward reaction rate constants, given in either
|
||||
* normal form or in exchange current density form.
|
||||
|
|
@ -467,48 +455,45 @@ protected:
|
|||
|
||||
//! Vector of chemical potentials for all species
|
||||
/*!
|
||||
* This vector contains a vector of chemical potentials for all of the species in the kinetics object
|
||||
* This vector contains a vector of chemical potentials for all of the
|
||||
* species in the kinetics object
|
||||
*
|
||||
* Length = m_kk. Units = J/kmol.
|
||||
*/
|
||||
vector_fp m_mu;
|
||||
|
||||
//! Vector of standard state electrochemical potentials modified by
|
||||
//! a standard concentration term.
|
||||
//! Vector of standard state electrochemical potentials modified by a
|
||||
//! standard concentration term.
|
||||
/*!
|
||||
* This vector contains a temporary vector of standard state electrochemical
|
||||
* potentials + RTln(Cs) for all of the species in the kinetics object
|
||||
*
|
||||
* In order to get the units correct for the concentration equilibrium
|
||||
* constant, each species needs to have an
|
||||
* RT ln(Cs) added to its contribution to the equilibrium constant
|
||||
* Cs is the standard concentration for the species. Frequently, for
|
||||
* solid species, Cs is equal to 1. However, for gases Cs is P/RT.
|
||||
*
|
||||
* Length = m_kk. Units = J/kmol.
|
||||
* constant, each species needs to have an RT ln(Cs) added to its
|
||||
* contribution to the equilibrium constant Cs is the standard concentration
|
||||
* for the species. Frequently, for solid species, Cs is equal to 1.
|
||||
* However, for gases Cs is P/RT. Length = m_kk. Units = J/kmol.
|
||||
*/
|
||||
vector_fp m_mu0_Kc;
|
||||
|
||||
//! Vector of phase electric potentials
|
||||
/*!
|
||||
* Temporary vector containing the potential of each phase in the kinetics
|
||||
* object.
|
||||
*
|
||||
* length = number of phases. Units = Volts.
|
||||
* 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.
|
||||
* 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;
|
||||
|
||||
//! Storage for the net electric energy change due to reaction.
|
||||
/*!
|
||||
* Length is number of reactions. It's used to store the
|
||||
* net electric potential energy change due to the reaction.
|
||||
* Length is number of reactions. It's used to store the net electric
|
||||
* potential energy change due to the reaction.
|
||||
*
|
||||
* deltaElectricEnergy_[jrxn] = sum_i ( F V_i z_i nu_ij)
|
||||
*/
|
||||
|
|
@ -516,8 +501,7 @@ protected:
|
|||
|
||||
//! Vector of raw activation energies for the reactions
|
||||
/*!
|
||||
* Units are in Kelvin.
|
||||
* Length is number of reactions.
|
||||
* Units are in Kelvin. Length is number of reactions.
|
||||
*/
|
||||
vector_fp m_E;
|
||||
|
||||
|
|
@ -527,15 +511,15 @@ protected:
|
|||
//! 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 object's surface
|
||||
* problem uncoupled from other surface phases.
|
||||
* be used to solve this single InterfaceKinetics object's surface problem
|
||||
* uncoupled from other surface phases.
|
||||
*/
|
||||
ImplicitSurfChem* m_integrator;
|
||||
|
||||
//! Electrochemical transfer coefficient for the forward direction
|
||||
/*!
|
||||
* Electrochemical transfer coefficient for all reactions that have transfer reactions
|
||||
* the reaction is given by m_ctrxn[i]
|
||||
* Electrochemical transfer coefficient for all reactions that have
|
||||
* transfer reactions the reaction is given by m_ctrxn[i]
|
||||
*/
|
||||
vector_fp m_beta;
|
||||
|
||||
|
|
@ -657,11 +641,11 @@ protected:
|
|||
|
||||
//! Vector of booleans indicating whether phases exist or not
|
||||
/*!
|
||||
* Vector of booleans indicating whether a phase exists or not. We use
|
||||
* this to set the ROP's so that unphysical things don't happen.
|
||||
* For example, a reaction can't go in the forwards direction if a
|
||||
* phase in which a reactant is present doesn't exist. Because InterfaceKinetics
|
||||
* deals with intrinsic quantities only normally, nowhere else is this extrinsic
|
||||
* Vector of booleans indicating whether a phase exists or not. We use this
|
||||
* to set the ROP's so that unphysical things don't happen. For example, a
|
||||
* reaction can't go in the forwards direction if a phase in which a
|
||||
* reactant is present doesn't exist. Because InterfaceKinetics deals with
|
||||
* intrinsic quantities only normally, nowhere else is this extrinsic
|
||||
* concept introduced except here.
|
||||
*
|
||||
* length = number of phases in the object. By default all phases exist.
|
||||
|
|
@ -678,8 +662,8 @@ protected:
|
|||
*/
|
||||
vector_int m_phaseIsStable;
|
||||
|
||||
//! Vector of vector of booleans indicating whether a phase participates in a
|
||||
//! reaction as a reactant
|
||||
//! Vector of vector of booleans indicating whether a phase participates in
|
||||
//! a reaction as a reactant
|
||||
/*!
|
||||
* m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p
|
||||
* participates in reaction j as a reactant.
|
||||
|
|
|
|||
|
|
@ -26,98 +26,84 @@ namespace Cantera
|
|||
/// @defgroup kineticsmgr Kinetics Managers
|
||||
/// @section kinmodman Models and Managers
|
||||
///
|
||||
/// A kinetics manager is a C++ class that implements a kinetics
|
||||
/// model; a kinetics model is a set of mathematical equation
|
||||
/// describing how various kinetic quantities are to be computed --
|
||||
/// reaction rates, species production rates, etc. Many different
|
||||
/// kinetics models might be defined to handle different types of
|
||||
/// kinetic processes. For example, one kinetics model might use
|
||||
/// expressions valid for elementary reactions in ideal gas
|
||||
/// mixtures. It might, for example, require the reaction orders
|
||||
/// to be integral and equal to the forward stoichiometric
|
||||
/// coefficients, require that each reaction be reversible with a
|
||||
/// reverse rate satisfying detailed balance, include
|
||||
/// pressure-dependent unimolecular reactions, etc. Another
|
||||
/// kinetics model might be designed for heterogeneous chemistry
|
||||
/// at interfaces, and might allow empirical reaction orders,
|
||||
/// coverage-dependent activation energies, irreversible
|
||||
/// reactions, and include effects of potential differences across
|
||||
/// A kinetics manager is a C++ class that implements a kinetics model; a
|
||||
/// kinetics model is a set of mathematical equation describing how various
|
||||
/// kinetic quantities are to be computed -- reaction rates, species production
|
||||
/// rates, etc. Many different kinetics models might be defined to handle
|
||||
/// different types of kinetic processes. For example, one kinetics model might
|
||||
/// use expressions valid for elementary reactions in ideal gas mixtures. It
|
||||
/// might, for example, require the reaction orders to be integral and equal to
|
||||
/// the forward stoichiometric coefficients, require that each reaction be
|
||||
/// reversible with a reverse rate satisfying detailed balance, include
|
||||
/// pressure-dependent unimolecular reactions, etc. Another kinetics model might
|
||||
/// be designed for heterogeneous chemistry at interfaces, and might allow
|
||||
/// empirical reaction orders, coverage-dependent activation energies,
|
||||
/// irreversible reactions, and include effects of potential differences across
|
||||
/// the interface on reaction rates.
|
||||
///
|
||||
/// A kinetics manager implements a kinetics model. Since the
|
||||
/// model equations may be complex and expensive to evaluate, a
|
||||
/// kinetics manager may adopt various strategies to 'manage' the
|
||||
/// computation and evaluate the expressions efficiently. For
|
||||
/// example, if there are rate coefficients or other quantities
|
||||
/// that depend only on temperature, a manager class may choose to
|
||||
/// store these quantities internally, and re-evaluate them only
|
||||
/// when the temperature has actually changed. Or a manager
|
||||
/// designed for use with reaction mechanisms with a few repeated
|
||||
/// activation energies might precompute the terms \f$ exp(-E/RT)
|
||||
/// \f$, instead of evaluating the exponential repeatedly for each
|
||||
/// reaction. There are many other possible 'management styles',
|
||||
/// each of which might be better suited to some reaction
|
||||
/// mechanisms than others.
|
||||
/// A kinetics manager implements a kinetics model. Since the model equations
|
||||
/// may be complex and expensive to evaluate, a kinetics manager may adopt
|
||||
/// various strategies to 'manage' the computation and evaluate the expressions
|
||||
/// efficiently. For example, if there are rate coefficients or other quantities
|
||||
/// that depend only on temperature, a manager class may choose to store these
|
||||
/// quantities internally, and re-evaluate them only when the temperature has
|
||||
/// actually changed. Or a manager designed for use with reaction mechanisms
|
||||
/// with a few repeated activation energies might precompute the terms \f$
|
||||
/// exp(-E/RT) \f$, instead of evaluating the exponential repeatedly for each
|
||||
/// reaction. There are many other possible 'management styles', each of which
|
||||
/// might be better suited to some reaction mechanisms than others.
|
||||
///
|
||||
/// But however a manager structures the internal computation, the
|
||||
/// tasks the manager class must perform are, for the most part,
|
||||
/// the same. It must be able to compute reaction rates, species
|
||||
/// production rates, equilibrium constants, etc. Therefore, all
|
||||
/// kinetics manager classes should have a common set of public
|
||||
/// methods, but differ in how they implement these methods.
|
||||
/// But however a manager structures the internal computation, the tasks the
|
||||
/// manager class must perform are, for the most part, the same. It must be able
|
||||
/// to compute reaction rates, species production rates, equilibrium constants,
|
||||
/// etc. Therefore, all kinetics manager classes should have a common set of
|
||||
/// public methods, but differ in how they implement these methods.
|
||||
///
|
||||
/// A kinetics manager computes reaction rates of progress,
|
||||
/// species production rates, equilibrium constants, and similar
|
||||
/// quantities for a reaction mechanism. All kinetics manager
|
||||
/// classes derive from class Kinetics, which defines a common
|
||||
/// public interface for all kinetics managers. Each derived class
|
||||
/// overloads the virtual methods of Kinetics to implement a
|
||||
/// particular kinetics model.
|
||||
/// A kinetics manager computes reaction rates of progress, species production
|
||||
/// rates, equilibrium constants, and similar quantities for a reaction
|
||||
/// mechanism. All kinetics manager classes derive from class Kinetics, which
|
||||
/// defines a common public interface for all kinetics managers. Each derived
|
||||
/// class overloads the virtual methods of Kinetics to implement a particular
|
||||
/// kinetics model.
|
||||
///
|
||||
/// For example, class GasKinetics implements reaction rate
|
||||
/// expressions appropriate for homogeneous reactions in ideal gas
|
||||
/// mixtures, and class InterfaceKinetics implements expressions
|
||||
/// appropriate for heterogeneous mechanisms at interfaces,
|
||||
/// including how to handle reactions involving charged species of
|
||||
/// phases with different electric potentials --- something that
|
||||
/// class GasKinetics doesn't deal with at all.
|
||||
/// For example, class GasKinetics implements reaction rate expressions
|
||||
/// appropriate for homogeneous reactions in ideal gas mixtures, and class
|
||||
/// InterfaceKinetics implements expressions appropriate for heterogeneous
|
||||
/// mechanisms at interfaces, including how to handle reactions involving
|
||||
/// charged species of phases with different electric potentials --- something
|
||||
/// that class GasKinetics doesn't deal with at all.
|
||||
///
|
||||
/// Many of the methods of class Kinetics write into arrays the
|
||||
/// values of some quantity for each species, for example the net
|
||||
/// production rate. These methods always write the results into
|
||||
/// flat arrays, ordered by phase in the order the phase was
|
||||
/// added, and within a phase in the order the species were added
|
||||
/// to the phase (which is the same ordering as in the input
|
||||
/// file). Example: suppose a heterogeneous mechanism involves
|
||||
/// three phases -- a bulk phase 'a', another bulk phase 'b', and
|
||||
/// the surface phase 'a:b' at the a/b interface. Phase 'a'
|
||||
/// contains 12 species, phase 'b' contains 3, and at the
|
||||
/// interface there are 5 adsorbed species defined in phase
|
||||
/// 'a:b'. Then methods like getNetProductionRates(doublereal* net)
|
||||
/// will write and output array of length 20, beginning at the location
|
||||
/// pointed to by 'net'. The first 12 values will be the net production
|
||||
/// rates for all 12 species of phase 'a' (even if some do not participate
|
||||
/// in the reactions), the next 3 will be for phase 'b', and finally the
|
||||
/// net production rates for the surface species will occupy the last
|
||||
/// 5 locations.
|
||||
/// Many of the methods of class Kinetics write into arrays the values of some
|
||||
/// quantity for each species, for example the net production rate. These
|
||||
/// methods always write the results into flat arrays, ordered by phase in the
|
||||
/// order the phase was added, and within a phase in the order the species were
|
||||
/// added to the phase (which is the same ordering as in the input file).
|
||||
/// Example: suppose a heterogeneous mechanism involves three phases -- a bulk
|
||||
/// phase 'a', another bulk phase 'b', and the surface phase 'a:b' at the a/b
|
||||
/// interface. Phase 'a' contains 12 species, phase 'b' contains 3, and at the
|
||||
/// interface there are 5 adsorbed species defined in phase 'a:b'. Then methods
|
||||
/// like getNetProductionRates(doublereal* net) will write and output array of
|
||||
/// length 20, beginning at the location pointed to by 'net'. The first 12
|
||||
/// values will be the net production rates for all 12 species of phase 'a'
|
||||
/// (even if some do not participate in the reactions), the next 3 will be for
|
||||
/// phase 'b', and finally the net production rates for the surface species will
|
||||
/// occupy the last 5 locations.
|
||||
/// @ingroup chemkinetics
|
||||
|
||||
|
||||
//! Public interface for kinetics managers.
|
||||
/*!
|
||||
* This class serves as a base class to derive 'kinetics
|
||||
* managers', which are classes that manage homogeneous chemistry
|
||||
* within one phase, or heterogeneous chemistry at one
|
||||
* interface. The virtual methods of this class are meant to be
|
||||
* overloaded in subclasses. The non-virtual methods perform
|
||||
* generic functions and are implemented in Kinetics. They should
|
||||
* not be overloaded. Only those methods required by a subclass
|
||||
* need to be overloaded; the rest will throw exceptions if
|
||||
* called.
|
||||
* This class serves as a base class to derive 'kinetics managers', which are
|
||||
* classes that manage homogeneous chemistry within one phase, or heterogeneous
|
||||
* chemistry at one interface. The virtual methods of this class are meant to be
|
||||
* overloaded in subclasses. The non-virtual methods perform generic functions
|
||||
* and are implemented in Kinetics. They should not be overloaded. Only those
|
||||
* methods required by a subclass need to be overloaded; the rest will throw
|
||||
* exceptions if called.
|
||||
*
|
||||
* When the nomenclature "kinetics species index" is used below,
|
||||
* this means that the species index ranges over all species in
|
||||
* all phases handled by the kinetics manager.
|
||||
* When the nomenclature "kinetics species index" is used below, this means that
|
||||
* the species index ranges over all species in all phases handled by the
|
||||
* kinetics manager.
|
||||
*
|
||||
* @ingroup kineticsmgr
|
||||
*/
|
||||
|
|
@ -164,8 +150,9 @@ public:
|
|||
* own the ThermoPhase objects. After a duplication, we need to point to
|
||||
* different ThermoPhase objects.
|
||||
*
|
||||
* We check that the ThermoPhase objects are aligned in the same order and have
|
||||
* the following identical properties to the ones that they are replacing:
|
||||
* We check that the ThermoPhase objects are aligned in the same order and
|
||||
* have the following identical properties to the ones that they are
|
||||
* replacing:
|
||||
*
|
||||
* - ThermoPhase::id()
|
||||
* - ThermoPhase::eosType()
|
||||
|
|
@ -178,9 +165,8 @@ public:
|
|||
|
||||
//! Identifies the kinetics manager type.
|
||||
/*!
|
||||
* Each class derived from Kinetics should overload this method to
|
||||
* return a unique integer. Standard values are defined in file
|
||||
* mix_defs.h.
|
||||
* Each class derived from Kinetics should overload this method to return
|
||||
* a unique integer. Standard values are defined in file mix_defs.h.
|
||||
*/
|
||||
virtual int type() const;
|
||||
|
||||
|
|
@ -212,10 +198,10 @@ public:
|
|||
//@{
|
||||
|
||||
/**
|
||||
* The number of phases participating in the reaction
|
||||
* mechanism. For a homogeneous reaction mechanism, this will
|
||||
* always return 1, but for a heterogeneous mechanism it will
|
||||
* return the total number of phases in the mechanism.
|
||||
* The number of phases participating in the reaction mechanism. For a
|
||||
* homogeneous reaction mechanism, this will always return 1, but for a
|
||||
* heterogeneous mechanism it will return the total number of phases in the
|
||||
* mechanism.
|
||||
*/
|
||||
size_t nPhases() const {
|
||||
return m_thermo.size();
|
||||
|
|
@ -231,13 +217,13 @@ public:
|
|||
void checkPhaseArraySize(size_t mm) const;
|
||||
|
||||
/**
|
||||
* Return the phase index of a phase in the list of phases
|
||||
* defined within the object.
|
||||
* Return the phase index of a phase in the list of phases defined within
|
||||
* the object.
|
||||
*
|
||||
* @param ph std::string name of the phase
|
||||
*
|
||||
* If a -1 is returned, then the phase is not defined in
|
||||
* the Kinetics object.
|
||||
* If a -1 is returned, then the phase is not defined in the Kinetics
|
||||
* object.
|
||||
*/
|
||||
size_t phaseIndex(const std::string& ph) {
|
||||
if (m_phaseindex.find(ph) == m_phaseindex.end()) {
|
||||
|
|
@ -261,7 +247,7 @@ public:
|
|||
* the phases in the list of phases represents the 2D interface or 1D edge
|
||||
* at which the reactions take place. This method returns the index of the
|
||||
* phase with the smallest spatial dimension (1, 2, or 3) among the list
|
||||
* of phases. If there is more than one, the index of the first one is
|
||||
* of phases. If there is more than one, the index of the first one is
|
||||
* returned. For homogeneous mechanisms, the value 0 is returned.
|
||||
*/
|
||||
size_t reactionPhaseIndex() {
|
||||
|
|
@ -270,7 +256,7 @@ public:
|
|||
|
||||
/**
|
||||
* This method returns a reference to the nth ThermoPhase object defined
|
||||
* in this kinetics mechanism. It is typically used so that member
|
||||
* in this kinetics mechanism. It is typically used so that member
|
||||
* functions of the ThermoPhase object may be called. For homogeneous
|
||||
* mechanisms, there is only one object, and this method can be called
|
||||
* without an argument to access it.
|
||||
|
|
@ -294,23 +280,19 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* The location of species k of phase n in species arrays.
|
||||
* Kinetics manager classes return species production rates in
|
||||
* flat arrays, with the species of each phases following one
|
||||
* another, in the order the phases were added. This method
|
||||
* is useful to find the value for a particular species of a
|
||||
* particular phase in arrays returned from methods like
|
||||
* getCreationRates that return an array of species-specific
|
||||
* quantities.
|
||||
* The location of species k of phase n in species arrays. Kinetics manager
|
||||
* classes return species production rates in flat arrays, with the species
|
||||
* of each phases following one another, in the order the phases were added.
|
||||
* This method is useful to find the value for a particular species of a
|
||||
* particular phase in arrays returned from methods like getCreationRates
|
||||
* that return an array of species-specific quantities.
|
||||
*
|
||||
* Example: suppose a heterogeneous mechanism involves three
|
||||
* phases. The first contains 12 species, the second 26, and
|
||||
* the third 3. Then species arrays must have size at least
|
||||
* 41, and positions 0 - 11 are the values for the species in
|
||||
* the first phase, positions 12 - 37 are the values for the
|
||||
* species in the second phase, etc. Then
|
||||
* kineticsSpeciesIndex(7, 0) = 7, kineticsSpeciesIndex(4, 1)
|
||||
* = 16, and kineticsSpeciesIndex(2, 2) = 40.
|
||||
* Example: suppose a heterogeneous mechanism involves three phases. The
|
||||
* first contains 12 species, the second 26, and the third 3. Then species
|
||||
* arrays must have size at least 41, and positions 0 - 11 are the values
|
||||
* for the species in the first phase, positions 12 - 37 are the values for
|
||||
* the species in the second phase, etc. Then kineticsSpeciesIndex(7, 0) =
|
||||
* 7, kineticsSpeciesIndex(4, 1) = 16, and kineticsSpeciesIndex(2, 2) = 40.
|
||||
*
|
||||
* @param k species index
|
||||
* @param n phase index for the species
|
||||
|
|
@ -323,7 +305,7 @@ public:
|
|||
/*!
|
||||
* k is an integer from 0 to ktot - 1, where ktot is the number of
|
||||
* species in the kinetics manager, which is the sum of the number of
|
||||
* species in all phases participating in the kinetics manager. If k is
|
||||
* species in all phases participating in the kinetics manager. If k is
|
||||
* out of bounds, the string "<unknown>" is returned.
|
||||
*
|
||||
* @param k species index
|
||||
|
|
@ -393,7 +375,7 @@ public:
|
|||
|
||||
//! Return the forward rates of progress of the reactions
|
||||
/*!
|
||||
* Forward rates of progress. Return the forward rates of
|
||||
* Forward rates of progress. Return the forward rates of
|
||||
* progress in array fwdROP, which must be dimensioned at
|
||||
* least as large as the total number of reactions.
|
||||
*
|
||||
|
|
@ -413,7 +395,7 @@ public:
|
|||
virtual void getRevRatesOfProgress(doublereal* revROP);
|
||||
|
||||
/**
|
||||
* Net rates of progress. Return the net (forward - reverse) rates of
|
||||
* Net rates of progress. Return the net (forward - reverse) rates of
|
||||
* progress in array netROP, which must be dimensioned at least as large
|
||||
* as the total number of reactions.
|
||||
*
|
||||
|
|
@ -439,17 +421,15 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* Change in species properties. Given an array of molar species
|
||||
* property values \f$ z_k, k = 1, \dots, K \f$, return the
|
||||
* array of reaction values
|
||||
* Change in species properties. Given an array of molar species property
|
||||
* values \f$ z_k, k = 1, \dots, K \f$, return the array of reaction values
|
||||
* \f[
|
||||
* \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I.
|
||||
* \f]
|
||||
* For example, if this method is called with the array of
|
||||
* standard-state molar Gibbs free energies for the species,
|
||||
* then the values returned in array \c deltaProperty would be
|
||||
* the standard-state Gibbs free energies of reaction for each
|
||||
* reaction.
|
||||
* For example, if this method is called with the array of standard-state
|
||||
* molar Gibbs free energies for the species, then the values returned in
|
||||
* array \c deltaProperty would be the standard-state Gibbs free energies of
|
||||
* reaction for each reaction.
|
||||
*
|
||||
* @param property Input vector of property value. Length: m_kk.
|
||||
* @param deltaProperty Output vector of deltaRxn. Length: nReactions().
|
||||
|
|
@ -512,7 +492,7 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* Return the vector of values for the reactions change in entropy. These
|
||||
* Return the vector of values for the reactions change in entropy. These
|
||||
* values depend upon the concentration of the solution.
|
||||
*
|
||||
* units = J kmol-1 Kelvin-1
|
||||
|
|
@ -525,9 +505,9 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* 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.
|
||||
* 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.
|
||||
*
|
||||
* units = J kmol-1
|
||||
*
|
||||
|
|
@ -539,9 +519,9 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* 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.
|
||||
* 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.
|
||||
*
|
||||
* units = J kmol-1
|
||||
*
|
||||
|
|
@ -553,9 +533,9 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* 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.
|
||||
* 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.
|
||||
*
|
||||
* units = J kmol-1 Kelvin-1
|
||||
*
|
||||
|
|
@ -580,19 +560,19 @@ public:
|
|||
virtual void getCreationRates(doublereal* cdot);
|
||||
|
||||
/**
|
||||
* Species destruction rates [kmol/m^3/s or 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. @see nTotalSpecies.
|
||||
* Species destruction rates [kmol/m^3/s or 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. @see nTotalSpecies.
|
||||
*
|
||||
* @param ddot Output vector of destruction rates. Length: m_kk.
|
||||
*/
|
||||
virtual void getDestructionRates(doublereal* ddot);
|
||||
|
||||
/**
|
||||
* Species net production rates [kmol/m^3/s or 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. @see nTotalSpecies.
|
||||
* Species net production rates [kmol/m^3/s or 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. @see nTotalSpecies.
|
||||
*
|
||||
* @param wdot Output vector of net production rates. Length: m_kk.
|
||||
*/
|
||||
|
|
@ -653,9 +633,8 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* Flag specifying the type of reaction. The legal values and
|
||||
* their meaning are specific to the particular kinetics
|
||||
* manager.
|
||||
* Flag specifying the type of reaction. The legal values and their meaning
|
||||
* are specific to the particular kinetics manager.
|
||||
*
|
||||
* @param i reaction index
|
||||
*/
|
||||
|
|
@ -664,9 +643,9 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* 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.
|
||||
* 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.
|
||||
*
|
||||
* @param i reaction index
|
||||
*/
|
||||
|
|
@ -809,7 +788,7 @@ public:
|
|||
/*!
|
||||
* These methods alter reaction rates. They are designed primarily for
|
||||
* carrying out sensitivity analysis, but may be used for any purpose
|
||||
* requiring dynamic alteration of rate constants. For each reaction, a
|
||||
* requiring dynamic alteration of rate constants. For each reaction, a
|
||||
* real-valued multiplier may be defined that multiplies the reaction rate
|
||||
* coefficient. The multiplier may be set to zero to completely remove a
|
||||
* reaction from the mechanism.
|
||||
|
|
@ -836,8 +815,8 @@ public:
|
|||
//@}
|
||||
|
||||
/**
|
||||
* Returns true if the kinetics manager has been properly
|
||||
* initialized and finalized.
|
||||
* Returns true if the kinetics manager has been properly initialized and
|
||||
* finalized.
|
||||
*/
|
||||
virtual bool ready() const {
|
||||
return false;
|
||||
|
|
@ -927,14 +906,12 @@ protected:
|
|||
//! m_thermo is a vector of pointers to ThermoPhase objects that are
|
||||
//! involved with this kinetics operator
|
||||
/*!
|
||||
* For homogeneous kinetics applications, this vector
|
||||
* will only have one entry. For interfacial reactions, this
|
||||
* vector will consist of multiple entries; some of them will
|
||||
* be surface phases, and the other ones will be bulk phases.
|
||||
* The order that the objects are listed determines the order
|
||||
* in which the species comprising each phase are listed in
|
||||
* the source term vector, originating from the reaction
|
||||
* mechanism.
|
||||
* For homogeneous kinetics applications, this vector will only have one
|
||||
* entry. For interfacial reactions, this vector will consist of multiple
|
||||
* entries; some of them will be surface phases, and the other ones will be
|
||||
* bulk phases. The order that the objects are listed determines the order
|
||||
* in which the species comprising each phase are listed in the source term
|
||||
* vector, originating from the reaction mechanism.
|
||||
*
|
||||
* Note that this kinetics object doesn't own these ThermoPhase objects
|
||||
* and is not responsible for creating or deleting them.
|
||||
|
|
@ -942,19 +919,16 @@ protected:
|
|||
std::vector<thermo_t*> m_thermo;
|
||||
|
||||
/**
|
||||
* m_start is a vector of integers specifying the beginning position
|
||||
* for the species vector for the n'th phase in the kinetics
|
||||
* class.
|
||||
* m_start is a vector of integers specifying the beginning position for the
|
||||
* species vector for the n'th phase in the kinetics class.
|
||||
*/
|
||||
std::vector<size_t> m_start;
|
||||
|
||||
/**
|
||||
* Mapping of the phase id, i.e., the id attribute in the XML
|
||||
* phase element to the position of the phase within the
|
||||
* kinetics object. Positions start with the value of 1. The
|
||||
* member function, phaseIndex() decrements by one before
|
||||
* returning the index value, so that missing phases return
|
||||
* -1.
|
||||
* Mapping of the phase id, i.e., the id attribute in the XML phase element
|
||||
* to the position of the phase within the kinetics object. Positions start
|
||||
* with the value of 1. The member function, phaseIndex() decrements by one
|
||||
* before returning the index value, so that missing phases return -1.
|
||||
*/
|
||||
std::map<std::string, size_t> m_phaseindex;
|
||||
|
||||
|
|
@ -963,8 +937,8 @@ protected:
|
|||
|
||||
//! Phase Index where reactions are assumed to be taking place
|
||||
/*!
|
||||
* We calculate this by assuming that the phase with the lowest
|
||||
* dimensionality is the phase where reactions are taking place.
|
||||
* We calculate this by assuming that the phase with the lowest
|
||||
* dimensionality is the phase where reactions are taking place.
|
||||
*/
|
||||
size_t m_rxnphase;
|
||||
|
||||
|
|
|
|||
|
|
@ -41,14 +41,12 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* Update the concentration-dependent parts of the rate
|
||||
* coefficient, if any. Used by class SurfaceArrhenius to
|
||||
* compute coverage-dependent * modifications to the Arrhenius
|
||||
* parameters. The array c should contain whatever data the
|
||||
* particular rate coefficient class needs to update its
|
||||
* rates. Note that this method does not return anything. To
|
||||
* get the updated rates, method update must be called after
|
||||
* the call to update_C.
|
||||
* Update the concentration-dependent parts of the rate coefficient, if any.
|
||||
* Used by class SurfaceArrhenius to compute coverage-dependent *
|
||||
* modifications to the Arrhenius parameters. The array c should contain
|
||||
* whatever data the particular rate coefficient class needs to update its
|
||||
* rates. Note that this method does not return anything. To get the
|
||||
* updated rates, method update must be called after the call to update_C.
|
||||
*/
|
||||
void update_C(const doublereal* c) {
|
||||
for (size_t i = 0; i != m_rates.size(); i++) {
|
||||
|
|
@ -57,12 +55,11 @@ public:
|
|||
}
|
||||
|
||||
/**
|
||||
* Write the rate coefficients into array values. Each
|
||||
* calculator writes one entry in values, at the location
|
||||
* specified by the reaction number when it was
|
||||
* installed. Note that nothing will be done for reactions
|
||||
* that have constant rates. The array values should be
|
||||
* preloaded with the constant rate coefficients.
|
||||
* Write the rate coefficients into array values. Each calculator writes one
|
||||
* entry in values, at the location specified by the reaction number when it
|
||||
* was installed. Note that nothing will be done for reactions that have
|
||||
* constant rates. The array values should be preloaded with the constant
|
||||
* rate coefficients.
|
||||
*/
|
||||
void update(doublereal T, doublereal logT, doublereal* values) {
|
||||
doublereal recipT = 1.0/T;
|
||||
|
|
|
|||
|
|
@ -85,8 +85,7 @@ public:
|
|||
typedef std::map<size_t, doublereal> rxn_path_map;
|
||||
|
||||
/**
|
||||
* Constructor. Construct a one-way path from
|
||||
* \c begin to \c end.
|
||||
* Constructor. Construct a one-way path from \c begin to \c end.
|
||||
*/
|
||||
Path(SpeciesNode* begin, SpeciesNode* end);
|
||||
|
||||
|
|
@ -194,12 +193,10 @@ public:
|
|||
/**
|
||||
* Export the reaction path diagram. This method writes to stream
|
||||
* \c s the commands for the 'dot' program in the \c GraphViz
|
||||
* package from AT&T. (GraphViz may be downloaded from
|
||||
* www.graphviz.org.)
|
||||
* package from AT&T. (GraphViz may be downloaded from www.graphviz.org.)
|
||||
*
|
||||
* To generate a postscript reaction path diagram from the
|
||||
* output of this method saved in file paths.dot, for example, give
|
||||
* the command:
|
||||
* To generate a postscript reaction path diagram from the output of this
|
||||
* method saved in file paths.dot, for example, give the command:
|
||||
* \code
|
||||
* dot -Tps paths.dot > paths.ps
|
||||
* \endcode
|
||||
|
|
@ -302,8 +299,7 @@ public:
|
|||
int build(Kinetics& s, const std::string& element, std::ostream& output,
|
||||
ReactionPathDiagram& r, bool quiet=false);
|
||||
|
||||
//! Analyze a reaction to determine which reactants lead to which
|
||||
//! products.
|
||||
//! Analyze a reaction to determine which reactants lead to which products.
|
||||
int findGroups(std::ostream& logfile, Kinetics& s);
|
||||
|
||||
void writeGroup(std::ostream& out, const Group& g);
|
||||
|
|
|
|||
|
|
@ -38,16 +38,16 @@ public:
|
|||
|
||||
/// Constructor.
|
||||
/// @param A pre-exponential. The unit system is
|
||||
/// (kmol, m, s). The actual units depend on the reaction
|
||||
/// order and the dimensionality (surface or bulk).
|
||||
/// (kmol, m, s). The actual units depend on the reaction
|
||||
/// order and the dimensionality (surface or bulk).
|
||||
/// @param b Temperature exponent. Non-dimensional.
|
||||
/// @param E Activation energy in temperature units. Kelvin.
|
||||
Arrhenius(doublereal A, doublereal b, doublereal E);
|
||||
|
||||
//! Update concentration-dependent parts of the rate coefficient.
|
||||
/*!
|
||||
* For this class, there are no
|
||||
* concentration-dependent parts, so this method does nothing.
|
||||
* For this class, there are no concentration-dependent parts, so this
|
||||
* method does nothing.
|
||||
*/
|
||||
void update_C(const doublereal* c) {
|
||||
}
|
||||
|
|
@ -62,9 +62,8 @@ public:
|
|||
/**
|
||||
* Update the value the rate constant.
|
||||
*
|
||||
* This function returns the actual value of the rate constant.
|
||||
* It can be safely called for negative values of the pre-exponential
|
||||
* factor.
|
||||
* This function returns the actual value of the rate constant. It can be
|
||||
* safely called for negative values of the pre-exponential factor.
|
||||
*/
|
||||
doublereal updateRC(doublereal logT, doublereal recipT) const {
|
||||
return m_A * std::exp(m_b*logT - m_E*recipT);
|
||||
|
|
|
|||
|
|
@ -55,10 +55,9 @@ namespace Cantera
|
|||
* this matrix for elementary reactions involving three or fewer product
|
||||
* molecules (or reactant molecules).
|
||||
*
|
||||
* To take advantage of this structure, reactions are divided into four
|
||||
* groups. These classes are designed to take advantage of this sparse
|
||||
* structure when computing quantities that can be written as matrix
|
||||
* multiplies.
|
||||
* To take advantage of this structure, reactions are divided into four groups.
|
||||
* These classes are designed to take advantage of this sparse structure when
|
||||
* computing quantities that can be written as matrix multiplies.
|
||||
*
|
||||
* They are designed to explicitly unroll loops over species or reactions for
|
||||
* Operations on reactions that require knowing the reaction stoichiometry.
|
||||
|
|
@ -118,11 +117,6 @@ namespace Cantera
|
|||
* by always assuming it is equal to one and then treating reactants and
|
||||
* products for a reaction separately. Bimolecular reactions involving the
|
||||
* identical species are treated as involving separate species.
|
||||
*
|
||||
* @internal This class should be upgraded to include cases where
|
||||
* real stoichiometric coefficients are used. Shouldn't be that
|
||||
* hard to do, and they occur in engineering simulations with some
|
||||
* regularity.
|
||||
*/
|
||||
|
||||
static doublereal ppow(doublereal x, doublereal order)
|
||||
|
|
@ -414,40 +408,40 @@ public:
|
|||
private:
|
||||
//! Length of the m_ic vector
|
||||
/*!
|
||||
* This is the number of species which participate in the reaction order
|
||||
* and stoichiometric coefficient vectors for the reactant or product description
|
||||
* of the reaction.
|
||||
* This is the number of species which participate in the reaction order
|
||||
* and stoichiometric coefficient vectors for the reactant or product
|
||||
* description of the reaction.
|
||||
*/
|
||||
size_t m_n;
|
||||
|
||||
//! ID of the reaction corresponding to this stoichiometric manager
|
||||
//! ID of the reaction corresponding to this stoichiometric manager
|
||||
/*!
|
||||
* This is used within the interface to select the array position to read and write to
|
||||
* Normally this is associated with the reaction number in an array of quantities indexed
|
||||
* by the reaction number, e.g., ROP[irxn].
|
||||
* This is used within the interface to select the array position to read
|
||||
* and write to Normally this is associated with the reaction number in an
|
||||
* array of quantities indexed by the reaction number, e.g., ROP[irxn].
|
||||
*/
|
||||
size_t m_rxn;
|
||||
|
||||
//! Vector of species which are involved with this stoichiometric manager calculations
|
||||
//! Vector of species which are involved with this stoichiometric manager
|
||||
//! calculations
|
||||
/*!
|
||||
* This is an integer list of species which participate in either the
|
||||
* reaction order matrix or the stoichiometric order matrix for this reaction, m_rxn.
|
||||
* reaction order matrix or the stoichiometric order matrix for this
|
||||
* reaction, m_rxn.
|
||||
*/
|
||||
std::vector<size_t> m_ic;
|
||||
|
||||
//! Reaction orders for the reaction
|
||||
/*!
|
||||
* This is either for the reactants or products.
|
||||
* Length = m_n
|
||||
* Species number, m_ic[n], has a reaction order of m_order[n].
|
||||
* This is either for the reactants or products. Length = m_n. Species
|
||||
* number, m_ic[n], has a reaction order of m_order[n].
|
||||
*/
|
||||
vector_fp m_order;
|
||||
|
||||
//! Stoichiometric coefficients for the reaction, reactant or product side.
|
||||
/*!
|
||||
* This is either for the reactants or products.
|
||||
* Length = m_n
|
||||
* Species number m_ic[m], has a stoichiometric coefficient of m_stoich[n].
|
||||
* This is either for the reactants or products. Length = m_n. Species
|
||||
* number m_ic[m], has a stoichiometric coefficient of m_stoich[n].
|
||||
*/
|
||||
vector_fp m_stoich;
|
||||
};
|
||||
|
|
@ -498,38 +492,36 @@ inline static void _decrementReactions(InputIter begin,
|
|||
}
|
||||
|
||||
/*
|
||||
* This class handles operations involving the stoichiometric
|
||||
* coefficients on one side of a reaction (reactant or product) for
|
||||
* a set of reactions comprising a reaction mechanism. This class is
|
||||
* used by class Kinetics, which contains three instances
|
||||
* of this class (one to handle operations on the reactions, one for
|
||||
* the products of reversible reactions, and one for the products of
|
||||
* irreversible reactions).
|
||||
* This class handles operations involving the stoichiometric coefficients on
|
||||
* one side of a reaction (reactant or product) for a set of reactions
|
||||
* comprising a reaction mechanism. This class is used by class Kinetics, which
|
||||
* contains three instances of this class (one to handle operations on the
|
||||
* reactions, one for the products of reversible reactions, and one for the
|
||||
* products of irreversible reactions).
|
||||
*
|
||||
* This class is designed for use with elementary reactions, or at
|
||||
* least ones with integral stoichiometric coefficients. Let \f$ M(i) \f$
|
||||
* be the number of molecules on the product or reactant side of
|
||||
* reaction number i.
|
||||
* This class is designed for use with elementary reactions, or at least ones
|
||||
* with integral stoichiometric coefficients. Let \f$ M(i) \f$ be the number of
|
||||
* molecules on the product or reactant side of reaction number i.
|
||||
* \f[
|
||||
* r_i = \sum_m^{M_i} s_{k_{m,i}}
|
||||
* \f]
|
||||
* To understand the operations performed by this class, let
|
||||
* \f$ N_{k,i}\f$ denote the stoichiometric coefficient of species k on
|
||||
* one side (reactant or product) in reaction i. Then \b N is a sparse
|
||||
* K by I matrix of stoichiometric coefficients.
|
||||
* To understand the operations performed by this class, let \f$ N_{k,i}\f$
|
||||
* denote the stoichiometric coefficient of species k on one side (reactant or
|
||||
* product) in reaction i. Then \b N is a sparse K by I matrix of stoichiometric
|
||||
* coefficients.
|
||||
*
|
||||
* The following matrix operations may be carried out with a vector
|
||||
* S of length K, and a vector R of length I:
|
||||
* The following matrix operations may be carried out with a vector S of length
|
||||
* K, and a vector R of length I:
|
||||
*
|
||||
* - \f$ S = S + N R\f$ (incrementSpecies)
|
||||
* - \f$ S = S - N R\f$ (decrementSpecies)
|
||||
* - \f$ R = R + N^T S \f$ (incrementReaction)
|
||||
* - \f$ R = R - N^T S \f$ (decrementReaction)
|
||||
*
|
||||
* The actual implementation, however, does not compute these
|
||||
* quantities by matrix multiplication. A faster algorithm is used
|
||||
* that makes use of the fact that the \b integer-valued N matrix is
|
||||
* very sparse, and the non-zero terms are small positive integers.
|
||||
* The actual implementation, however, does not compute these quantities by
|
||||
* matrix multiplication. A faster algorithm is used that makes use of the fact
|
||||
* that the \b integer-valued N matrix is very sparse, and the non-zero terms
|
||||
* are small positive integers.
|
||||
* \f[
|
||||
* S_k = R_{i1} + \dots + R_{iM}
|
||||
* \f]
|
||||
|
|
@ -543,26 +535,25 @@ public:
|
|||
/**
|
||||
* Constructor for the StoichManagerN class.
|
||||
*
|
||||
* @internal Consider adding defaulted entries here that supply
|
||||
* the total number of reactions in the mechanism and the total
|
||||
* number of species in the species list. Then, we could use those
|
||||
* numbers to provide error checks during the construction of the
|
||||
* object. Those numbers would also provide some clarity to the
|
||||
* purpose and utility of this class.
|
||||
* @internal Consider adding defaulted entries here that supply the total
|
||||
* number of reactions in the mechanism and the total number of species
|
||||
* in the species list. Then, we could use those numbers to provide
|
||||
* error checks during the construction of the object. Those numbers
|
||||
* would also provide some clarity to the purpose and utility of this
|
||||
* class.
|
||||
*
|
||||
* DGG - the problem is that the number of reactions and species
|
||||
* are not known initially.
|
||||
* DGG - the problem is that the number of reactions and species are not
|
||||
* known initially.
|
||||
*/
|
||||
StoichManagerN() {
|
||||
}
|
||||
|
||||
/**
|
||||
* Add a single reaction to the list of reactions that this
|
||||
* stoichiometric manager object handles.
|
||||
* Add a single reaction to the list of reactions that this stoichiometric
|
||||
* manager object handles.
|
||||
*
|
||||
* This function is the same as the add() function below. However,
|
||||
* the order of each species in the power list expression is
|
||||
* set to one automatically.
|
||||
* This function is the same as the add() function below. However, the order
|
||||
* of each species in the power list expression is set to one automatically.
|
||||
*/
|
||||
void add(size_t rxn, const std::vector<size_t>& k) {
|
||||
vector_fp order(k.size(), 1.0);
|
||||
|
|
@ -578,19 +569,18 @@ public:
|
|||
//! Add a single reaction to the list of reactions that this
|
||||
//! stoichiometric manager object handles.
|
||||
/*!
|
||||
* @param rxn Reaction index of the current reaction. This is used
|
||||
* as an index into vectors which have length n_total_rxn.
|
||||
* @param k This is a vector of integer values specifying the
|
||||
* species indices. The length of this vector species
|
||||
* the number of different species in the description.
|
||||
* The value of the entries are the species indices.
|
||||
* These are used as indexes into vectors which have
|
||||
* length n_total_species.
|
||||
* @param order This is a vector of the same length as vector k.
|
||||
* The order is used for the routine power(), which produces
|
||||
* a power law expression involving the species vector.
|
||||
* @param stoich This is used to handle fractional stoichiometric coefficients
|
||||
* on the product side of irreversible reactions.
|
||||
* @param rxn Reaction index of the current reaction. This is used as an
|
||||
* index into vectors which have length n_total_rxn.
|
||||
* @param k This is a vector of integer values specifying the species
|
||||
* indices. The length of this vector species the number of different
|
||||
* species in the description. The value of the entries are the species
|
||||
* indices. These are used as indexes into vectors which have length
|
||||
* n_total_species.
|
||||
* @param order This is a vector of the same length as vector k. The order
|
||||
* is used for the routine power(), which produces a power law
|
||||
* expression involving the species vector.
|
||||
* @param stoich This is used to handle fractional stoichiometric
|
||||
* coefficients on the product side of irreversible reactions.
|
||||
*/
|
||||
void add(size_t rxn, const std::vector<size_t>& k, const vector_fp& order,
|
||||
const vector_fp& stoich) {
|
||||
|
|
|
|||
|
|
@ -21,10 +21,9 @@ namespace Cantera
|
|||
|
||||
//! Install information about reactions into the kinetics object, kin.
|
||||
/*!
|
||||
* At this point, parent usually refers to the phase XML element.
|
||||
* One of the children of this element is reactionArray,
|
||||
* the element which determines where in the XML file to
|
||||
* look up the reaction rate data.
|
||||
* At this point, parent usually refers to the phase XML element. One of the
|
||||
* children of this element is reactionArray, the element which determines
|
||||
* where in the XML file to look up the reaction rate data.
|
||||
*
|
||||
* @param p parent XML phase element
|
||||
* @param kin Kinetics object to install reactions into
|
||||
|
|
|
|||
|
|
@ -50,8 +50,8 @@ const int FALLOFF_RXN = 4;
|
|||
const int PLOG_RXN = 5;
|
||||
|
||||
/**
|
||||
* A general gas-phase pressure-dependent reaction where k(T,P) is defined in terms of
|
||||
* a bivariate Chebyshev polynomial.
|
||||
* A general gas-phase pressure-dependent reaction where k(T,P) is defined in
|
||||
* terms of a bivariate Chebyshev polynomial.
|
||||
*/
|
||||
const int CHEBYSHEV_RXN = 6;
|
||||
|
||||
|
|
|
|||
|
|
@ -70,9 +70,8 @@ namespace Cantera
|
|||
* 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.
|
||||
* 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
|
||||
|
|
@ -204,9 +203,8 @@ private:
|
|||
//! 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.
|
||||
* 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.
|
||||
|
|
@ -308,10 +306,10 @@ private:
|
|||
const doublereal* CSolnSPOld, const bool do_time,
|
||||
const doublereal deltaT);
|
||||
|
||||
//! Pointer to the manager of the implicit surface chemistry problem
|
||||
//! 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.
|
||||
* This object actually calls the current object. Thus, we are providing a
|
||||
* loop-back functionality here.
|
||||
*/
|
||||
ImplicitSurfChem* m_SurfChemPtr;
|
||||
|
||||
|
|
@ -363,7 +361,7 @@ private:
|
|||
*/
|
||||
std::vector<size_t> m_nSpeciesSurfPhase;
|
||||
|
||||
//! Vector of surface phase pointers
|
||||
//! Vector of surface phase pointers
|
||||
/*!
|
||||
* This is created during the constructor
|
||||
* Length is equal to the number of surface phases, m_numSurfPhases
|
||||
|
|
@ -480,10 +478,10 @@ private:
|
|||
//! 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)
|
||||
//! Weights for the residual norm calculation. length MAX(1, m_neq)
|
||||
vector_fp m_wtResid;
|
||||
|
||||
//! Weights for the species concentrations norm calculation
|
||||
//! Weights for the species concentrations norm calculation
|
||||
/*!
|
||||
* length MAX(1, m_neq)
|
||||
*/
|
||||
|
|
@ -491,12 +489,11 @@ private:
|
|||
|
||||
//! 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.
|
||||
* 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)
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -1,8 +1,7 @@
|
|||
/**
|
||||
* @file AqueousKinetics.cpp
|
||||
*
|
||||
* Homogeneous kinetics in an aqueous phase, either condensed
|
||||
* or dilute in salts
|
||||
* Homogeneous kinetics in an aqueous phase, either condensed or dilute in salts
|
||||
*/
|
||||
/*
|
||||
* Copyright (2006) Sandia Corporation. Under the terms of
|
||||
|
|
|
|||
|
|
@ -86,11 +86,8 @@ void BulkKinetics::getDeltaSSEntropy(doublereal* deltaS)
|
|||
|
||||
void BulkKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible)
|
||||
{
|
||||
/*
|
||||
* go get the forward rate constants. -> note, we don't
|
||||
* really care about speed or redundancy in these
|
||||
* informational routines.
|
||||
*/
|
||||
// go get the forward rate constants. -> note, we don't really care about
|
||||
// speed or redundancy in these informational routines.
|
||||
getFwdRateConstants(krev);
|
||||
|
||||
if (doIrreversible) {
|
||||
|
|
|
|||
|
|
@ -1,7 +1,5 @@
|
|||
/**
|
||||
* @file GasKinetics.cpp
|
||||
*
|
||||
* Homogeneous kinetics in ideal gases
|
||||
* @file GasKinetics.cpp Homogeneous kinetics in ideal gases
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
|
@ -271,8 +269,8 @@ bool GasKinetics::addReaction(shared_ptr<Reaction> r)
|
|||
|
||||
void GasKinetics::addFalloffReaction(FalloffReaction& r)
|
||||
{
|
||||
// install high and low rate coeff calculators
|
||||
// and extend the high and low rate coeff value vectors
|
||||
// install high and low rate coeff calculators and extend the high and low
|
||||
// rate coeff value vectors
|
||||
m_falloff_high_rates.install(m_nfall, r.high_rate);
|
||||
m_rfn_high.push_back(0.0);
|
||||
m_falloff_low_rates.install(m_nfall, r.low_rate);
|
||||
|
|
|
|||
|
|
@ -1,13 +1,10 @@
|
|||
/**
|
||||
* @file Group.cpp
|
||||
*
|
||||
* Implementation file for the Group class used in reaction path analysis.
|
||||
* @file Group.cpp Implementation file for the Group class used in reaction path
|
||||
* analysis.
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
||||
// reaction path analysis support
|
||||
|
||||
#include "cantera/kinetics/Group.h"
|
||||
#include <iostream>
|
||||
|
||||
|
|
|
|||
|
|
@ -170,20 +170,13 @@ void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
|
|||
doublereal timeScaleOverride)
|
||||
{
|
||||
int ifunc;
|
||||
/*
|
||||
* set bulkFunc
|
||||
* -> We assume that the bulk concentrations are constant.
|
||||
*/
|
||||
// set bulkFunc. We assume that the bulk concentrations are constant.
|
||||
int bulkFunc = BULK_ETCH;
|
||||
/*
|
||||
* time scale - time over which to integrate equations
|
||||
*/
|
||||
// time scale - time over which to integrate equations
|
||||
doublereal time_scale = timeScaleOverride;
|
||||
if (!m_surfSolver) {
|
||||
m_surfSolver.reset(new solveSP(this, bulkFunc));
|
||||
/*
|
||||
* set ifunc, which sets the algorithm.
|
||||
*/
|
||||
// set ifunc, which sets the algorithm.
|
||||
ifunc = SFLUX_INITIALIZE;
|
||||
} else {
|
||||
ifunc = SFLUX_RESIDUAL;
|
||||
|
|
@ -194,24 +187,20 @@ void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
|
|||
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
|
||||
*/
|
||||
// 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(m_concSpecies.data());
|
||||
InterfaceKinetics* ik = m_vecKinPtrs[0];
|
||||
ThermoPhase& tp = ik->thermo(0);
|
||||
doublereal TKelvin = tp.temperature();
|
||||
doublereal 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.
|
||||
*/
|
||||
|
||||
// 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);
|
||||
}
|
||||
|
|
@ -219,11 +208,8 @@ void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
|
|||
doublereal reltol = 1.0E-6;
|
||||
doublereal 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.
|
||||
*/
|
||||
// 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 (size_t k = 0; k < m_nv; k++) {
|
||||
if (m_concSpecies[k] < 0.0) {
|
||||
|
|
|
|||
|
|
@ -47,17 +47,13 @@ InterfaceKinetics::~InterfaceKinetics()
|
|||
|
||||
InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right)
|
||||
{
|
||||
/*
|
||||
* Call the assignment operator
|
||||
*/
|
||||
// Call the assignment operator
|
||||
operator=(right);
|
||||
}
|
||||
|
||||
InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
|
||||
{
|
||||
/*
|
||||
* Check for self assignment.
|
||||
*/
|
||||
// Check for self assignment.
|
||||
if (this == &right) {
|
||||
return *this;
|
||||
}
|
||||
|
|
@ -147,8 +143,8 @@ void InterfaceKinetics::_update_rates_T()
|
|||
m_rates.update(T, m_logtemp, m_rfn.data());
|
||||
applyStickingCorrection(m_rfn.data());
|
||||
|
||||
// If we need to do conversions between exchange current density formulation and regular formulation
|
||||
// (either way) do it here.
|
||||
// If we need to do conversions between exchange current density
|
||||
// formulation and regular formulation (either way) do it here.
|
||||
if (m_has_exchange_current_density_formulation) {
|
||||
convertExchangeCurrentDensityFormulation(m_rfn.data());
|
||||
}
|
||||
|
|
@ -173,18 +169,16 @@ void InterfaceKinetics::_update_rates_phi()
|
|||
}
|
||||
}
|
||||
|
||||
// Updates the internal variables m_actConc and m_conc
|
||||
void InterfaceKinetics::_update_rates_C()
|
||||
{
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
const ThermoPhase* tp = m_thermo[n];
|
||||
/*
|
||||
* We call the getActivityConcentrations function of each
|
||||
* ThermoPhase class that makes up this kinetics object to
|
||||
* obtain the generalized concentrations for species within that
|
||||
* class. This is collected in the vector m_conc. m_start[]
|
||||
* are integer indices for that vector denoting the start of the
|
||||
* species for each phase.
|
||||
* We call the getActivityConcentrations function of each ThermoPhase
|
||||
* class that makes up this kinetics object to obtain the generalized
|
||||
* concentrations for species within that class. This is collected in
|
||||
* the vector m_conc. m_start[] are integer indices for that vector
|
||||
* denoting the start of the species for each phase.
|
||||
*/
|
||||
tp->getActivityConcentrations(m_actConc.data() + m_start[n]);
|
||||
|
||||
|
|
@ -206,8 +200,9 @@ void InterfaceKinetics::updateKc()
|
|||
|
||||
if (m_nrev > 0) {
|
||||
/*
|
||||
* Get the vector of standard state electrochemical potentials for species in the Interfacial
|
||||
* kinetics object and store it in m_mu0[] and m_mu0_Kc[]
|
||||
* Get the vector of standard state electrochemical potentials for
|
||||
* species in the Interfacial kinetics object and store it in m_mu0[]
|
||||
* and m_mu0_Kc[]
|
||||
*/
|
||||
updateMu0();
|
||||
doublereal rrt = 1.0 / (GasConstant * thermo(0).temperature());
|
||||
|
|
@ -235,10 +230,6 @@ void InterfaceKinetics::updateMu0()
|
|||
_update_rates_phi();
|
||||
|
||||
updateExchangeCurrentQuantities();
|
||||
/*
|
||||
* Get the vector of standard state electrochemical potentials for species in the Interfacial
|
||||
* kinetics object and store it in m_mu0[] and in m_mu0_Kc[]
|
||||
*/
|
||||
size_t nsp, ik = 0;
|
||||
size_t np = nPhases();
|
||||
for (size_t n = 0; n < np; n++) {
|
||||
|
|
@ -300,20 +291,16 @@ void InterfaceKinetics::getEquilibriumConstants(doublereal* kc)
|
|||
|
||||
void InterfaceKinetics::updateExchangeCurrentQuantities()
|
||||
{
|
||||
/*
|
||||
* Calculate:
|
||||
* - m_StandardConc[]
|
||||
* - m_ProdStandConcReac[]
|
||||
* - m_deltaG0[]
|
||||
* - m_mu0[]
|
||||
*/
|
||||
// Calculate:
|
||||
// - m_StandardConc[]
|
||||
// - m_ProdStanConcReac[]
|
||||
// - m_deltaG0[]
|
||||
// - m_mu0[]
|
||||
|
||||
/*
|
||||
* First collect vectors of the standard Gibbs free energies of the
|
||||
* species and the standard concentrations
|
||||
* - m_mu0
|
||||
* - m_StandardConc
|
||||
*/
|
||||
// First collect vectors of the standard Gibbs free energies of the
|
||||
// species and the standard concentrations
|
||||
// - m_mu0
|
||||
// - m_StandardConc
|
||||
size_t ik = 0;
|
||||
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
|
|
@ -346,27 +333,23 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf)
|
|||
}
|
||||
}
|
||||
|
||||
// Compute the change in electrical potential energy for each
|
||||
// reaction. This will only be non-zero if a potential
|
||||
// difference is present.
|
||||
// Compute the change in electrical potential energy for each reaction. This
|
||||
// will only be non-zero if a potential difference is present.
|
||||
getReactionDelta(m_pot.data(), deltaElectricEnergy_.data());
|
||||
|
||||
// Modify the reaction rates. Only modify those with a
|
||||
// non-zero activation energy. Below we decrease the
|
||||
// activation energy below zero but in some debug modes
|
||||
// we print out a warning message about this.
|
||||
/*
|
||||
* NOTE, there is some discussion about this point.
|
||||
* Should we decrease the activation energy below zero?
|
||||
* I don't think this has been decided in any definitive way.
|
||||
* The treatment below is numerically more stable, however.
|
||||
*/
|
||||
// Modify the reaction rates. Only modify those with a non-zero activation
|
||||
// energy. Below we decrease the activation energy below zero but in some
|
||||
// debug modes we print out a warning message about this.
|
||||
|
||||
// NOTE, there is some discussion about this point. Should we decrease the
|
||||
// activation energy below zero? I don't think this has been decided in any
|
||||
// definitive way. The treatment below is numerically more stable, however.
|
||||
doublereal eamod;
|
||||
for (size_t i = 0; i < m_beta.size(); i++) {
|
||||
size_t irxn = m_ctrxn[i];
|
||||
|
||||
// If we calculate the BV form directly, we don't add the voltage correction to the
|
||||
// forward reaction rate constants.
|
||||
// If we calculate the BV form directly, we don't add the voltage
|
||||
// correction to the forward reaction rate constants.
|
||||
if (m_ctrxn_BVform[i] == 0) {
|
||||
eamod = m_beta[i] * deltaElectricEnergy_[irxn];
|
||||
if (eamod != 0.0) {
|
||||
|
|
@ -379,18 +362,21 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf)
|
|||
void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* const kfwd)
|
||||
{
|
||||
updateExchangeCurrentQuantities();
|
||||
// Loop over all reactions which are defined to have a voltage transfer coefficient that
|
||||
// affects the activity energy for the reaction
|
||||
// Loop over all reactions which are defined to have a voltage transfer
|
||||
// coefficient that affects the activity energy for the reaction
|
||||
for (size_t i = 0; i < m_ctrxn.size(); i++) {
|
||||
size_t irxn = m_ctrxn[i];
|
||||
|
||||
// Determine whether the reaction rate constant is in an exchange current density formulation format.
|
||||
// Determine whether the reaction rate constant is in an exchange
|
||||
// current density formulation format.
|
||||
int iECDFormulation = m_ctrxn_ecdf[i];
|
||||
if (iECDFormulation) {
|
||||
// If the BV form is to be converted into the normal form then we go through this process.
|
||||
// If it isn't to be converted, then we don't go through this process.
|
||||
// If the BV form is to be converted into the normal form then we go
|
||||
// through this process. If it isn't to be converted, then we don't
|
||||
// go through this process.
|
||||
//
|
||||
// We need to have the straight chemical reaction rate constant to come out of this calculation.
|
||||
// We need to have the straight chemical reaction rate constant to
|
||||
// come out of this calculation.
|
||||
if (m_ctrxn_BVform[i] == 0) {
|
||||
// Calculate the term and modify the forward reaction
|
||||
double tmp = exp(- m_beta[i] * m_deltaG0[irxn] / thermo(0).RT());
|
||||
|
|
@ -402,12 +388,13 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con
|
|||
} else {
|
||||
// kfwd[] is the chemical reaction rate constant
|
||||
//
|
||||
// If we are to calculate the BV form directly, then we will do the reverse.
|
||||
// We will calculate the exchange current density formulation here and
|
||||
// substitute it.
|
||||
// If we are to calculate the BV form directly, then we will do the
|
||||
// reverse. We will calculate the exchange current density
|
||||
// formulation here and substitute it.
|
||||
if (m_ctrxn_BVform[i] != 0) {
|
||||
// Calculate the term and modify the forward reaction rate constant so that
|
||||
// it's in the exchange current density formulation format
|
||||
// Calculate the term and modify the forward reaction rate
|
||||
// constant so that it's in the exchange current density
|
||||
// formulation format
|
||||
double tmp = exp(m_beta[i] * m_deltaG0[irxn] * thermo(0).RT());
|
||||
double tmp2 = m_ProdStanConcReac[irxn];
|
||||
tmp *= Faraday * tmp2;
|
||||
|
|
@ -443,7 +430,8 @@ void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversibl
|
|||
|
||||
void InterfaceKinetics::updateROP()
|
||||
{
|
||||
// evaluate rate constants and equilibrium constants at temperature and phi (electric potential)
|
||||
// evaluate rate constants and equilibrium constants at temperature and phi
|
||||
// (electric potential)
|
||||
_update_rates_T();
|
||||
// get updated activities (rates updated below)
|
||||
_update_rates_C();
|
||||
|
|
@ -470,15 +458,17 @@ void InterfaceKinetics::updateROP()
|
|||
// the forward rates of progress.
|
||||
m_reactantStoich.multiply(m_actConc.data(), m_ropf.data());
|
||||
|
||||
// For reversible reactions, multiply ropr by the activity concentration products
|
||||
// For reversible reactions, multiply ropr by the activity concentration
|
||||
// products
|
||||
m_revProductStoich.multiply(m_actConc.data(), m_ropr.data());
|
||||
|
||||
// Fix up these calculations for cases where the above formalism doesn't hold
|
||||
// Fix up these calculations for cases where the above formalism doesn't hold
|
||||
double OCV = 0.0;
|
||||
for (size_t jrxn = 0; jrxn != nReactions(); ++jrxn) {
|
||||
if (reactionType(jrxn) == BUTLERVOLMER_RXN) {
|
||||
// OK, the reaction rate constant contains the current density rate constant calculation
|
||||
// the rxnstoich calculation contained the dependence of the current density on the activity concentrations
|
||||
// OK, the reaction rate constant contains the current density rate
|
||||
// constant calculation the rxnstoich calculation contained the
|
||||
// dependence of the current density on the activity concentrations
|
||||
// We finish up with the ROP calculation
|
||||
//
|
||||
// Calculate the overpotential of the reaction
|
||||
|
|
@ -494,11 +484,10 @@ void InterfaceKinetics::updateROP()
|
|||
m_ropnet[j] = m_ropf[j] - m_ropr[j];
|
||||
}
|
||||
|
||||
/*
|
||||
* For reactions involving multiple phases, we must check that the phase
|
||||
* being consumed actually exists. This is particularly important for
|
||||
* phases that are stoichiometric phases containing one species with a unity activity
|
||||
*/
|
||||
// For reactions involving multiple phases, we must check that the phase
|
||||
// being consumed actually exists. This is particularly important for phases
|
||||
// that are stoichiometric phases containing one species with a unity
|
||||
// activity
|
||||
if (m_phaseExistsCheck) {
|
||||
for (size_t j = 0; j != nReactions(); ++j) {
|
||||
if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) {
|
||||
|
|
@ -547,10 +536,8 @@ void InterfaceKinetics::updateROP()
|
|||
|
||||
void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
|
||||
{
|
||||
/*
|
||||
* Get the chemical potentials of the species in the all of the phases used in the
|
||||
* kinetics mechanism
|
||||
*/
|
||||
// Get the chemical potentials of the species in the all of the phases used
|
||||
// in the kinetics mechanism
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
m_thermo[n]->getChemPotentials(m_mu.data() + m_start[n]);
|
||||
}
|
||||
|
|
@ -566,107 +553,82 @@ void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
|
|||
|
||||
void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM)
|
||||
{
|
||||
/*
|
||||
* Get the chemical potentials of the species
|
||||
*/
|
||||
// Get the chemical potentials of the species
|
||||
size_t np = nPhases();
|
||||
for (size_t n = 0; n < np; n++) {
|
||||
thermo(n).getElectrochemPotentials(m_grt.data() + m_start[n]);
|
||||
}
|
||||
/*
|
||||
* Use the stoichiometric manager to find deltaG for each
|
||||
* reaction.
|
||||
*/
|
||||
|
||||
// Use the stoichiometric manager to find deltaG for each reaction.
|
||||
getReactionDelta(m_grt.data(), deltaM);
|
||||
}
|
||||
|
||||
void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH)
|
||||
{
|
||||
/*
|
||||
* Get the partial molar enthalpy of all species
|
||||
*/
|
||||
// Get the partial molar enthalpy of all species
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
thermo(n).getPartialMolarEnthalpies(m_grt.data() + m_start[n]);
|
||||
}
|
||||
/*
|
||||
* Use the stoichiometric manager to find deltaG for each
|
||||
* reaction.
|
||||
*/
|
||||
|
||||
// Use the stoichiometric manager to find deltaH for each reaction.
|
||||
getReactionDelta(m_grt.data(), deltaH);
|
||||
}
|
||||
|
||||
void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS)
|
||||
{
|
||||
/*
|
||||
* Get the partial molar entropy of all species in all of
|
||||
* the phases
|
||||
*/
|
||||
// Get the partial molar entropy of all species in all of the phases
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
thermo(n).getPartialMolarEntropies(m_grt.data() + m_start[n]);
|
||||
}
|
||||
/*
|
||||
* Use the stoichiometric manager to find deltaS for each
|
||||
* reaction.
|
||||
*/
|
||||
|
||||
// Use the stoichiometric manager to find deltaS for each reaction.
|
||||
getReactionDelta(m_grt.data(), deltaS);
|
||||
}
|
||||
|
||||
void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS)
|
||||
{
|
||||
/*
|
||||
* Get the standard state chemical potentials of the species.
|
||||
* This is the array of chemical potentials at unit activity
|
||||
* We define these here as the chemical potentials of the pure
|
||||
* species at the temperature and pressure of the solution.
|
||||
*/
|
||||
// Get the standard state chemical potentials of the species. This is the
|
||||
// array of chemical potentials at unit activity We define these here as the
|
||||
// chemical potentials of the pure species at the temperature and pressure
|
||||
// of the solution.
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
thermo(n).getStandardChemPotentials(m_mu0.data() + m_start[n]);
|
||||
}
|
||||
/*
|
||||
* Use the stoichiometric manager to find deltaG for each
|
||||
* reaction.
|
||||
*/
|
||||
|
||||
// Use the stoichiometric manager to find deltaG for each reaction.
|
||||
getReactionDelta(m_mu0.data(), deltaGSS);
|
||||
}
|
||||
|
||||
void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
|
||||
{
|
||||
/*
|
||||
* Get the standard state enthalpies of the species.
|
||||
* This is the array of chemical potentials at unit activity
|
||||
* We define these here as the enthalpies of the pure
|
||||
* species at the temperature and pressure of the solution.
|
||||
*/
|
||||
// Get the standard state enthalpies of the species. This is the array of
|
||||
// chemical potentials at unit activity We define these here as the
|
||||
// enthalpies of the pure species at the temperature and pressure of the
|
||||
// solution.
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
thermo(n).getEnthalpy_RT(m_grt.data() + m_start[n]);
|
||||
}
|
||||
for (size_t k = 0; k < m_kk; k++) {
|
||||
m_grt[k] *= thermo(0).RT();
|
||||
}
|
||||
/*
|
||||
* Use the stoichiometric manager to find deltaG for each
|
||||
* reaction.
|
||||
*/
|
||||
|
||||
// Use the stoichiometric manager to find deltaH for each reaction.
|
||||
getReactionDelta(m_grt.data(), deltaH);
|
||||
}
|
||||
|
||||
void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
|
||||
{
|
||||
/*
|
||||
* Get the standard state entropy of the species.
|
||||
* We define these here as the entropies of the pure
|
||||
* species at the temperature and pressure of the solution.
|
||||
*/
|
||||
// Get the standard state entropy of the species. We define these here as
|
||||
// the entropies of the pure species at the temperature and pressure of the
|
||||
// solution.
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
thermo(n).getEntropy_R(m_grt.data() + m_start[n]);
|
||||
}
|
||||
for (size_t k = 0; k < m_kk; k++) {
|
||||
m_grt[k] *= GasConstant;
|
||||
}
|
||||
/*
|
||||
* Use the stoichiometric manager to find deltaS for each
|
||||
* reaction.
|
||||
*/
|
||||
|
||||
// Use the stoichiometric manager to find deltaS for each reaction.
|
||||
getReactionDelta(m_grt.data(), deltaS);
|
||||
}
|
||||
|
||||
|
|
@ -945,9 +907,7 @@ void InterfaceKinetics::solvePseudoSteadyStateProblem(
|
|||
m_integrator->initialize();
|
||||
}
|
||||
m_integrator->setIOFlag(m_ioFlag);
|
||||
/*
|
||||
* New direct method to go here
|
||||
*/
|
||||
// New direct method to go here
|
||||
m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
|
||||
}
|
||||
|
||||
|
|
@ -1062,9 +1022,10 @@ void InterfaceKinetics::applyStickingCorrection(double* kf)
|
|||
|
||||
void EdgeKinetics::finalize()
|
||||
{
|
||||
// Note we can't call the Interface::finalize() routine because we need to check for a dimension of 1 below.
|
||||
// Therefore, we have to malloc room in arrays that would normally be
|
||||
// handled by the InterfaceKinetics::finalize() call.
|
||||
// Note we can't call the Interface::finalize() routine because we need to
|
||||
// check for a dimension of 1 below. Therefore, we have to malloc room in
|
||||
// arrays that would normally be handled by the
|
||||
// InterfaceKinetics::finalize() call.
|
||||
Kinetics::finalize();
|
||||
|
||||
size_t safe_reaction_size = std::max<size_t>(nReactions(), 1);
|
||||
|
|
|
|||
|
|
@ -30,17 +30,13 @@ Kinetics::~Kinetics() {}
|
|||
|
||||
Kinetics::Kinetics(const Kinetics& right)
|
||||
{
|
||||
/*
|
||||
* Call the assignment operator
|
||||
*/
|
||||
// Call the assignment operator
|
||||
*this = right;
|
||||
}
|
||||
|
||||
Kinetics& Kinetics::operator=(const Kinetics& right)
|
||||
{
|
||||
/*
|
||||
* Check for self assignment.
|
||||
*/
|
||||
// Check for self assignment.
|
||||
if (this == &right) {
|
||||
return *this;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -22,17 +22,15 @@ std::mutex KineticsFactory::kinetics_mutex;
|
|||
Kinetics* KineticsFactory::newKinetics(XML_Node& phaseData,
|
||||
vector<ThermoPhase*> th)
|
||||
{
|
||||
/*
|
||||
* Look for a child of the XML element phase called
|
||||
* "kinetics". It has an attribute name "model".
|
||||
* Store the value of that attribute in the variable kintype
|
||||
*/
|
||||
// Look for a child of the XML element phase called "kinetics". It has an
|
||||
// attribute name "model". Store the value of that attribute in the variable
|
||||
// kintype
|
||||
string kintype = phaseData.child("kinetics")["model"];
|
||||
|
||||
// Create a kinetics object of the desired type
|
||||
Kinetics* k = newKinetics(kintype);
|
||||
// Now that we have the kinetics manager, we can
|
||||
// import the reaction mechanism into it.
|
||||
// Now that we have the kinetics manager, we can import the reaction
|
||||
// mechanism into it.
|
||||
importKinetics(phaseData, th, k);
|
||||
|
||||
// Return the pointer to the kinetics manager
|
||||
|
|
|
|||
|
|
@ -226,8 +226,8 @@ void ReactionPathDiagram::exportToDot(ostream& s)
|
|||
|
||||
// draw paths representing net flows
|
||||
if (flow_type == NetFlow) {
|
||||
// if no scale was specified, normalize
|
||||
// net flows by the maximum net flow
|
||||
// if no scale was specified, normalize net flows by the maximum net
|
||||
// flow
|
||||
if (scale <= 0.0) {
|
||||
for (i1 = 0; i1 < nNodes(); i1++) {
|
||||
k1 = m_speciesNumber[i1];
|
||||
|
|
@ -256,8 +256,8 @@ void ReactionPathDiagram::exportToDot(ostream& s)
|
|||
flx = 0.0;
|
||||
}
|
||||
if (flx != 0.0) {
|
||||
// set beginning and end of the path based on the
|
||||
// sign of the net flow
|
||||
// set beginning and end of the path based on the sign of
|
||||
// the net flow
|
||||
if (flx > 0.0) {
|
||||
kbegin = k1;
|
||||
kend = k2;
|
||||
|
|
@ -268,8 +268,8 @@ void ReactionPathDiagram::exportToDot(ostream& s)
|
|||
flxratio = -flx/flmax;
|
||||
}
|
||||
|
||||
// write out path specification if the net flow
|
||||
// is greater than the threshold
|
||||
// write out path specification if the net flow is greater
|
||||
// than the threshold
|
||||
if (flxratio >= threshold) {
|
||||
// make nodes visible
|
||||
node(kbegin)->visible = true;
|
||||
|
|
@ -445,8 +445,7 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s)
|
|||
size_t kp0 = m_prod[i][0];
|
||||
size_t kp1 = m_prod[i][1];
|
||||
|
||||
// references to the Group objects representing the
|
||||
// reactants
|
||||
// references to the Group objects representing the reactants
|
||||
const Group& r0 = m_sgroup[kr0];
|
||||
const Group& r1 = m_sgroup[kr1];
|
||||
const Group& p0 = m_sgroup[kp0];
|
||||
|
|
@ -557,10 +556,10 @@ void ReactionPathBuilder::findElements(Kinetics& kin)
|
|||
for (size_t m = 0; m < nel; m++) {
|
||||
ename = p->elementName(m);
|
||||
|
||||
// if no entry is found for this element name, then
|
||||
// it is a new element. In this case, add the name
|
||||
// to the list of names, increment the element count,
|
||||
// and add an entry to the name->(index+1) map.
|
||||
// if no entry is found for this element name, then it is a new
|
||||
// element. In this case, add the name to the list of names,
|
||||
// increment the element count, and add an entry to the
|
||||
// name->(index+1) map.
|
||||
if (m_enamemap.find(ename) == m_enamemap.end()) {
|
||||
m_enamemap[ename] = m_nel + 1;
|
||||
m_elementSymbols.push_back(ename);
|
||||
|
|
@ -597,8 +596,8 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin)
|
|||
m_ns = kin.nTotalSpecies();
|
||||
m_nr = kin.nReactions();
|
||||
|
||||
// all reactants / products, even ones appearing on both sides
|
||||
// of the reaction
|
||||
// all reactants / products, even ones appearing on both sides of the
|
||||
// reaction
|
||||
vector<vector<size_t> > allProducts(m_nr);
|
||||
vector<vector<size_t> > allReactants(m_nr);
|
||||
for (size_t i = 0; i < m_nr; i++) {
|
||||
|
|
@ -614,7 +613,6 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin)
|
|||
|
||||
// m_reac and m_prod exclude indices for species that appear on
|
||||
// both sides of the reaction, so that the diagram contains no loops.
|
||||
|
||||
m_reac.resize(m_nr);
|
||||
m_prod.resize(m_nr);
|
||||
m_ropf.resize(m_nr);
|
||||
|
|
@ -628,8 +626,8 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin)
|
|||
map<size_t, int> net;
|
||||
|
||||
for (size_t i = 0; i < m_nr; i++) {
|
||||
// construct the lists of reactant and product indices, not
|
||||
// including molecules that appear on both sides.
|
||||
// construct the lists of reactant and product indices, not including
|
||||
// molecules that appear on both sides.
|
||||
m_reac[i].clear();
|
||||
m_prod[i].clear();
|
||||
net.clear();
|
||||
|
|
@ -658,10 +656,9 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin)
|
|||
|
||||
size_t nrnet = m_reac[i].size();
|
||||
|
||||
// compute number of atoms of each element in each reaction,
|
||||
// excluding molecules that appear on both sides of the
|
||||
// reaction. We only need to compute this for the reactants,
|
||||
// since the elements are conserved.
|
||||
// compute number of atoms of each element in each reaction, excluding
|
||||
// molecules that appear on both sides of the reaction. We only need to
|
||||
// compute this for the reactants, since the elements are conserved.
|
||||
for (n = 0; n < nrnet; n++) {
|
||||
k = m_reac[i][n];
|
||||
for (size_t m = 0; m < m_nel; m++) {
|
||||
|
|
@ -680,13 +677,12 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin)
|
|||
m_sgroup[j] = Group(comp);
|
||||
}
|
||||
|
||||
// determine whether or not the reaction is "determinate", meaning
|
||||
// that there is no ambiguity about which reactant is the source for
|
||||
// any element in any product. This is false if more than one
|
||||
// reactant contains a given element, *and* more than one product
|
||||
// contains the element. In this case, additional information is
|
||||
// needed to determine the partitioning of the reactant atoms of
|
||||
// that element among the products.
|
||||
// determine whether or not the reaction is "determinate", meaning that
|
||||
// there is no ambiguity about which reactant is the source for any element
|
||||
// in any product. This is false if more than one reactant contains a given
|
||||
// element, *and* more than one product contains the element. In this case,
|
||||
// additional information is needed to determine the partitioning of the
|
||||
// reactant atoms of that element among the products.
|
||||
int nar, nap;
|
||||
for (size_t i = 0; i < m_nr; i++) {
|
||||
nr = m_reac[i].size();
|
||||
|
|
@ -790,19 +786,17 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element,
|
|||
revlabel += " (+ M)";
|
||||
}
|
||||
|
||||
// calculate the flow only for pairs that are
|
||||
// not the same species, both contain atoms of
|
||||
// element m, and both are allowed to appear in
|
||||
// the diagram
|
||||
// calculate the flow only for pairs that are not the same
|
||||
// species, both contain atoms of element m, and both are
|
||||
// allowed to appear in the diagram
|
||||
if ((kkr != kkp) && (m_atoms(kkr,m) > 0
|
||||
&& m_atoms(kkp,m) > 0)
|
||||
&& status[kkr] >= 0 && status[kkp] >= 0) {
|
||||
// if neither species contains the full
|
||||
// number of atoms of element m in the
|
||||
// reaction, then we must consider the
|
||||
// type of reaction to determine which
|
||||
// reactant species was the source of a
|
||||
// given m-atom in the product
|
||||
// if neither species contains the full number of atoms
|
||||
// of element m in the reaction, then we must consider
|
||||
// the type of reaction to determine which reactant
|
||||
// species was the source of a given m-atom in the
|
||||
// product
|
||||
if ((m_atoms(kkp,m) < m_elatoms(m, i)) &&
|
||||
(m_atoms(kkr,m) < m_elatoms(m, i))) {
|
||||
map<size_t, map<size_t, Group> >& g = m_transfer[i];
|
||||
|
|
@ -825,12 +819,11 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element,
|
|||
}
|
||||
}
|
||||
} else {
|
||||
// no ambiguity about where the m-atoms come
|
||||
// from or go to. Either all reactant m atoms
|
||||
// end up in one product, or only one reactant
|
||||
// contains all the m-atoms. In either case,
|
||||
// the number of atoms transferred is given by
|
||||
// the same expression.
|
||||
// no ambiguity about where the m-atoms come from or
|
||||
// go to. Either all reactant m atoms end up in one
|
||||
// product, or only one reactant contains all the
|
||||
// m-atoms. In either case, the number of atoms
|
||||
// transferred is given by the same expression.
|
||||
f = m_atoms(kkp,m) * m_atoms(kkr,m) / m_elatoms(m, i);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -28,46 +28,38 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin,
|
|||
std::string default_phase, bool check_for_duplicates)
|
||||
{
|
||||
int itot = 0;
|
||||
/*
|
||||
* Search the children of the phase element for the
|
||||
* XML element named reactionArray. If we can't find it,
|
||||
* then return signaling having not found any reactions.
|
||||
* Apparently, we allow multiple reactionArray elements here
|
||||
* Each one will be processed sequentially, with the
|
||||
* end result being purely additive.
|
||||
*/
|
||||
|
||||
// Search the children of the phase element for the XML element named
|
||||
// reactionArray. If we can't find it, then return signaling having not
|
||||
// found any reactions. Apparently, we allow multiple reactionArray elements
|
||||
// here Each one will be processed sequentially, with the end result being
|
||||
// purely additive.
|
||||
vector<XML_Node*> rarrays = p.getChildren("reactionArray");
|
||||
if (rarrays.empty()) {
|
||||
kin.finalize();
|
||||
return false;
|
||||
}
|
||||
for (size_t n = 0; n < rarrays.size(); n++) {
|
||||
/*
|
||||
* Go get a reference to the current XML element,
|
||||
* reactionArray. We will process this element now.
|
||||
*/
|
||||
// Go get a reference to the current XML element, reactionArray. We will
|
||||
// process this element now.
|
||||
const XML_Node& rxns = *rarrays[n];
|
||||
/*
|
||||
* The reactionArray element has an attribute called,
|
||||
* datasrc. The value of the attribute is the XML
|
||||
* element comprising the top of the
|
||||
* tree of reactions for the phase.
|
||||
* Find this datasrc element starting with the root
|
||||
* of the current XML node.
|
||||
*/
|
||||
|
||||
// The reactionArray element has an attribute called, datasrc. The value
|
||||
// of the attribute is the XML element comprising the top of the tree of
|
||||
// reactions for the phase. Find this datasrc element starting with the
|
||||
// root of the current XML node.
|
||||
const XML_Node* rdata = get_XML_Node(rxns["datasrc"], &rxns.root());
|
||||
/*
|
||||
* If the reactionArray element has a child element named "skip", and
|
||||
* if the attribute of skip called "species" has a value of "undeclared",
|
||||
* we will set rxnrule.skipUndeclaredSpecies to 'true'. rxnrule is
|
||||
* passed to the routine that parses each individual reaction so that
|
||||
* the parser will skip all reactions containing an undefined species
|
||||
* without throwing an error.
|
||||
*
|
||||
* Similarly, an attribute named "third_bodies" with the value of
|
||||
* "undeclared" will skip undeclared third body efficiencies (while
|
||||
* retaining the reaction and any other efficiencies).
|
||||
*/
|
||||
|
||||
// If the reactionArray element has a child element named "skip", and if
|
||||
// the attribute of skip called "species" has a value of "undeclared",
|
||||
// we will set rxnrule.skipUndeclaredSpecies to 'true'. rxnrule is
|
||||
// passed to the routine that parses each individual reaction so that
|
||||
// the parser will skip all reactions containing an undefined species
|
||||
// without throwing an error.
|
||||
//
|
||||
// Similarly, an attribute named "third_bodies" with the value of
|
||||
// "undeclared" will skip undeclared third body efficiencies (while
|
||||
// retaining the reaction and any other efficiencies).
|
||||
if (rxns.hasChild("skip")) {
|
||||
const XML_Node& sk = rxns.child("skip");
|
||||
if (sk["species"] == "undeclared") {
|
||||
|
|
@ -77,11 +69,10 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin,
|
|||
kin.skipUndeclaredThirdBodies(true);
|
||||
}
|
||||
}
|
||||
/*
|
||||
* Search for child elements called include. We only include
|
||||
* a reaction if it's tagged by one of the include fields.
|
||||
* Or, we include all reactions if there are no include fields.
|
||||
*/
|
||||
|
||||
// Search for child elements called include. We only include a reaction
|
||||
// if it's tagged by one of the include fields. Or, we include all
|
||||
// reactions if there are no include fields.
|
||||
vector<XML_Node*> incl = rxns.getChildren("include");
|
||||
vector<XML_Node*> allrxns = rdata->getChildren("reaction");
|
||||
// if no 'include' directive, then include all reactions
|
||||
|
|
@ -113,11 +104,10 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin,
|
|||
if (iwild != string::npos) {
|
||||
rxid = rxid.substr(0,iwild);
|
||||
}
|
||||
/*
|
||||
* To decide whether the reaction is included or not
|
||||
* we do a lexical min max and operation. This
|
||||
* sometimes has surprising results.
|
||||
*/
|
||||
|
||||
// To decide whether the reaction is included or not we
|
||||
// do a lexical min max and operation. This sometimes
|
||||
// has surprising results.
|
||||
if ((rxid >= imin) && (rxid <= imax)) {
|
||||
kin.addReaction(newReaction(*r));
|
||||
++itot;
|
||||
|
|
@ -131,10 +121,9 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin,
|
|||
if (check_for_duplicates) {
|
||||
kin.checkDuplicates();
|
||||
}
|
||||
/*
|
||||
* Finalize the installation of the kinetics, now that we know
|
||||
* the true number of reactions in the mechanism, itot.
|
||||
*/
|
||||
|
||||
// Finalize the installation of the kinetics, now that we know the true
|
||||
// number of reactions in the mechanism, itot.
|
||||
kin.finalize();
|
||||
return true;
|
||||
}
|
||||
|
|
@ -159,10 +148,9 @@ bool importKinetics(const XML_Node& phase, std::vector<ThermoPhase*> th,
|
|||
}
|
||||
}
|
||||
|
||||
// if other phases are involved in the reaction mechanism,
|
||||
// they must be listed in a 'phaseArray' child
|
||||
// element. Homogeneous mechanisms do not need to include a
|
||||
// phaseArray element.
|
||||
// if other phases are involved in the reaction mechanism, they must be
|
||||
// listed in a 'phaseArray' child element. Homogeneous mechanisms do not
|
||||
// need to include a phaseArray element.
|
||||
vector<string> phase_ids;
|
||||
if (phase.hasChild("phaseArray")) {
|
||||
const XML_Node& pa = phase.child("phaseArray");
|
||||
|
|
@ -200,9 +188,8 @@ bool importKinetics(const XML_Node& phase, std::vector<ThermoPhase*> th,
|
|||
}
|
||||
}
|
||||
|
||||
// allocates arrays, etc. Must be called after the phases have
|
||||
// been added to 'kin', so that the number of species in each
|
||||
// phase is known.
|
||||
// allocates arrays, etc. Must be called after the phases have been added to
|
||||
// 'kin', so that the number of species in each phase is known.
|
||||
k->init();
|
||||
|
||||
// Install the reactions.
|
||||
|
|
@ -218,23 +205,18 @@ bool buildSolutionFromXML(XML_Node& root, const std::string& id,
|
|||
return false;
|
||||
}
|
||||
|
||||
/*
|
||||
* Fill in the ThermoPhase object by querying the
|
||||
* const XML_Node tree located at x.
|
||||
*/
|
||||
// Fill in the ThermoPhase object by querying the const XML_Node tree
|
||||
// located at x.
|
||||
importPhase(*x, th);
|
||||
/*
|
||||
* Create a vector of ThermoPhase pointers of length 1
|
||||
* having the current th ThermoPhase as the entry.
|
||||
*/
|
||||
|
||||
// Create a vector of ThermoPhase pointers of length 1 having the current th
|
||||
// ThermoPhase as the entry.
|
||||
std::vector<ThermoPhase*> phases(1);
|
||||
phases[0] = th;
|
||||
/*
|
||||
* Fill in the kinetics object k, by querying the
|
||||
* const XML_Node tree located by x. The source terms and
|
||||
* eventually the source term vector will be constructed
|
||||
* from the list of ThermoPhases in the vector, phases.
|
||||
*/
|
||||
|
||||
// Fill in the kinetics object k, by querying the const XML_Node tree
|
||||
// located by x. The source terms and eventually the source term vector will
|
||||
// be constructed from the list of ThermoPhases in the vector, phases.
|
||||
importKinetics(*x, phases, kin);
|
||||
return true;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -16,16 +16,12 @@ using namespace std;
|
|||
namespace Cantera
|
||||
{
|
||||
|
||||
/***************************************************************************
|
||||
* STATIC ROUTINES DEFINED IN THIS FILE
|
||||
***************************************************************************/
|
||||
// STATIC ROUTINES DEFINED IN THIS FILE
|
||||
|
||||
static doublereal calc_damping(doublereal* x, doublereal* dx, size_t dim, int*);
|
||||
static doublereal calcWeightedNorm(const doublereal [], const doublereal dx[], size_t);
|
||||
|
||||
/***************************************************************************
|
||||
* solveSP Class Definitions
|
||||
***************************************************************************/
|
||||
// solveSP Class Definitions
|
||||
|
||||
solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) :
|
||||
m_SurfChemPtr(surfChemPtr),
|
||||
|
|
@ -68,9 +64,7 @@ solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) :
|
|||
m_nSpeciesSurfPhase.push_back(nsp);
|
||||
m_numTotSurfSpecies += nsp;
|
||||
}
|
||||
/*
|
||||
* We rely on ordering to figure things out
|
||||
*/
|
||||
// We rely on ordering to figure things out
|
||||
m_numBulkPhasesSS = 0;
|
||||
|
||||
if (bulkFunc == BULK_DEPOSITION) {
|
||||
|
|
@ -130,33 +124,29 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
EXTRA_ACCURACY *= 0.001;
|
||||
}
|
||||
int info = 0;
|
||||
int label_t=-1; /* Species IDs for time control */
|
||||
int label_d = -1; /* Species IDs for damping control */
|
||||
int label_t=-1; // Species IDs for time control
|
||||
int label_d = -1; // Species IDs for damping control
|
||||
int label_t_old=-1;
|
||||
doublereal label_factor = 1.0;
|
||||
int iter=0; // iteration number on numlinear solver
|
||||
int iter_max=1000; // maximum number of nonlinear iterations
|
||||
doublereal deltaT = 1.0E-10; // Delta time step
|
||||
doublereal damp=1.0, tmp;
|
||||
// Weighted L2 norm of the residual. Currently, this is only
|
||||
// used for IO purposes. It doesn't control convergence.
|
||||
// Weighted L2 norm of the residual. Currently, this is only used for IO
|
||||
// purposes. It doesn't control convergence.
|
||||
doublereal resid_norm;
|
||||
doublereal inv_t = 0.0;
|
||||
doublereal t_real = 0.0, update_norm = 1.0E6;
|
||||
bool do_time = false, not_converged = true;
|
||||
m_ioflag = std::min(m_ioflag, 1);
|
||||
|
||||
/*
|
||||
* Set the initial value of the do_time parameter
|
||||
*/
|
||||
// Set the initial value of the do_time parameter
|
||||
if (ifunc == SFLUX_INITIALIZE || ifunc == SFLUX_TRANSIENT) {
|
||||
do_time = true;
|
||||
}
|
||||
|
||||
/*
|
||||
* Store the initial guess for the surface problem in the soln vector,
|
||||
* CSoln, and in an separate vector CSolnInit.
|
||||
*/
|
||||
// Store the initial guess for the surface problem in the soln vector,
|
||||
// CSoln, and in an separate vector CSolnInit.
|
||||
size_t loc = 0;
|
||||
for (size_t n = 0; n < m_numSurfPhases; n++) {
|
||||
SurfPhase* sf_ptr = m_ptrsSurfPhase[n];
|
||||
|
|
@ -177,39 +167,28 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
print_header(m_ioflag, ifunc, time_scale, true, reltol, abstol);
|
||||
}
|
||||
|
||||
/*
|
||||
* Quick return when there isn't a surface problem to solve
|
||||
*/
|
||||
// Quick return when there isn't a surface problem to solve
|
||||
if (m_neq == 0) {
|
||||
not_converged = false;
|
||||
update_norm = 0.0;
|
||||
}
|
||||
|
||||
/* ------------------------------------------------------------------
|
||||
* Start of Newton's method
|
||||
* ------------------------------------------------------------------
|
||||
*/
|
||||
// Start of Newton's method
|
||||
while (not_converged && iter < iter_max) {
|
||||
iter++;
|
||||
/*
|
||||
* Store previous iteration's solution in the old solution vector
|
||||
*/
|
||||
// Store previous iteration's solution in the old solution vector
|
||||
std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPOld.begin());
|
||||
|
||||
/*
|
||||
* Evaluate the largest surface species for each surface phase every
|
||||
* 5 iterations.
|
||||
*/
|
||||
// Evaluate the largest surface species for each surface phase every
|
||||
// 5 iterations.
|
||||
if (iter%5 == 4) {
|
||||
evalSurfLarge(m_CSolnSP.data());
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate the value of the time step
|
||||
* - heuristics to stop large oscillations in deltaT
|
||||
*/
|
||||
// Calculate the value of the time step
|
||||
// - heuristics to stop large oscillations in deltaT
|
||||
if (do_time) {
|
||||
/* don't hurry increase in time step at the same time as damping */
|
||||
// don't hurry increase in time step at the same time as damping
|
||||
if (damp < 1.0) {
|
||||
label_factor = 1.0;
|
||||
}
|
||||
|
|
@ -224,9 +203,7 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
inv_t = tmp;
|
||||
}
|
||||
|
||||
/*
|
||||
* Check end condition
|
||||
*/
|
||||
// Check end condition
|
||||
if (ifunc == SFLUX_TRANSIENT) {
|
||||
tmp = t_real + 1.0/inv_t;
|
||||
if (tmp > time_scale) {
|
||||
|
|
@ -234,44 +211,34 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
}
|
||||
}
|
||||
} else {
|
||||
/* make steady state calc a step of 1 million seconds to
|
||||
prevent singular Jacobians for some pathological cases */
|
||||
// make steady state calc a step of 1 million seconds to prevent
|
||||
// singular Jacobians for some pathological cases
|
||||
inv_t = 1.0e-6;
|
||||
}
|
||||
deltaT = 1.0/inv_t;
|
||||
|
||||
/*
|
||||
* Call the routine to numerically evaluation the Jacobian
|
||||
* and residual for the current iteration.
|
||||
*/
|
||||
// Call the routine to numerically evaluation the Jacobian and residual
|
||||
// for the current iteration.
|
||||
resjac_eval(m_Jac, m_resid.data(), m_CSolnSP.data(),
|
||||
m_CSolnSPOld.data(), do_time, deltaT);
|
||||
|
||||
/*
|
||||
* Calculate the weights. Make sure the calculation is carried
|
||||
* out on the first iteration.
|
||||
*/
|
||||
// Calculate the weights. Make sure the calculation is carried out on
|
||||
// the first iteration.
|
||||
if (iter%4 == 1) {
|
||||
calcWeights(m_wtSpecies.data(), m_wtResid.data(),
|
||||
m_Jac, m_CSolnSP.data(), abstol, reltol);
|
||||
}
|
||||
|
||||
/*
|
||||
* Find the weighted norm of the residual
|
||||
*/
|
||||
// Find the weighted norm of the residual
|
||||
resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
|
||||
|
||||
/*
|
||||
* Solve Linear system. The solution is in resid[]
|
||||
*/
|
||||
// Solve Linear system. The solution is in resid[]
|
||||
info = m_Jac.factor();
|
||||
if (info==0) {
|
||||
m_Jac.solve(&m_resid[0]);
|
||||
} else {
|
||||
/*
|
||||
* Force convergence if residual is small to avoid
|
||||
* "nan" results from the linear solve.
|
||||
*/
|
||||
// Force convergence if residual is small to avoid "nan" results
|
||||
// from the linear solve.
|
||||
if (m_ioflag) {
|
||||
writelogf("solveSurfSS: Zero pivot, assuming converged: %g (%d)\n",
|
||||
resid_norm, info);
|
||||
|
|
@ -280,7 +247,7 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
m_resid[jcol] = 0.0;
|
||||
}
|
||||
|
||||
/* print out some helpful info */
|
||||
// print out some helpful info
|
||||
if (m_ioflag > 1) {
|
||||
writelog("-----\n");
|
||||
writelogf("solveSurfProb: iter %d t_real %g delta_t %g\n\n",
|
||||
|
|
@ -299,24 +266,17 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculate the Damping factor needed to keep all unknowns
|
||||
* between 0 and 1, and not allow too large a change (factor of 2)
|
||||
* in any unknown.
|
||||
*/
|
||||
// Calculate the Damping factor needed to keep all unknowns between 0
|
||||
// and 1, and not allow too large a change (factor of 2) in any unknown.
|
||||
damp = calc_damping(m_CSolnSP.data(), m_resid.data(), m_neq, &label_d);
|
||||
|
||||
/*
|
||||
* Calculate the weighted norm of the update vector
|
||||
* Here, resid is the delta of the solution, in concentration
|
||||
* units.
|
||||
*/
|
||||
// Calculate the weighted norm of the update vector Here, resid is the
|
||||
// delta of the solution, in concentration units.
|
||||
update_norm = calcWeightedNorm(m_wtSpecies.data(),
|
||||
m_resid.data(), m_neq);
|
||||
/*
|
||||
* Update the solution vector and real time
|
||||
* Crop the concentrations to zero.
|
||||
*/
|
||||
|
||||
// Update the solution vector and real time Crop the concentrations to
|
||||
// zero.
|
||||
for (size_t irow = 0; irow < m_neq; irow++) {
|
||||
m_CSolnSP[irow] -= damp * m_resid[irow];
|
||||
}
|
||||
|
|
@ -348,12 +308,10 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
(resid_norm > EXTRA_ACCURACY));
|
||||
}
|
||||
}
|
||||
} /* End of Newton's Method while statement */
|
||||
} // End of Newton's Method while statement
|
||||
|
||||
/*
|
||||
* End Newton's method. If not converged, print error message and
|
||||
* recalculate sdot's at equal site fractions.
|
||||
*/
|
||||
// End Newton's method. If not converged, print error message and
|
||||
// recalculate sdot's at equal site fractions.
|
||||
if (not_converged && m_ioflag) {
|
||||
writelog("#$#$#$# Error in solveSP $#$#$#$ \n");
|
||||
writelogf("Newton iter on surface species did not converge, "
|
||||
|
|
@ -361,11 +319,8 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
writelog("Continuing anyway\n");
|
||||
}
|
||||
|
||||
/*
|
||||
* Decide on what to return in the solution vector
|
||||
* - right now, will always return the last solution
|
||||
* no matter how bad
|
||||
*/
|
||||
// Decide on what to return in the solution vector. Right now, will always
|
||||
// return the last solution no matter how bad
|
||||
if (m_ioflag) {
|
||||
fun_eval(m_resid.data(), m_CSolnSP.data(), m_CSolnSPOld.data(),
|
||||
false, deltaT);
|
||||
|
|
@ -374,9 +329,7 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
|
|||
update_norm, resid_norm, do_time, true);
|
||||
}
|
||||
|
||||
/*
|
||||
* Return with the appropriate flag
|
||||
*/
|
||||
// Return with the appropriate flag
|
||||
if (update_norm > 1.0) {
|
||||
return -1;
|
||||
}
|
||||
|
|
@ -437,18 +390,15 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln,
|
|||
doublereal sd = 0.0;
|
||||
doublereal grRate;
|
||||
if (m_numSurfPhases > 0) {
|
||||
/*
|
||||
* update the surface concentrations with the input surface
|
||||
* concentration vector
|
||||
*/
|
||||
// update the surface concentrations with the input surface
|
||||
// concentration vector
|
||||
updateState(CSoln);
|
||||
/*
|
||||
* Get the net production rates of all of the species in the
|
||||
* surface kinetics mechanism
|
||||
*
|
||||
* HKM Should do it here for all kinetics objects so that
|
||||
* bulk will eventually work.
|
||||
*/
|
||||
|
||||
// Get the net production rates of all of the species in the
|
||||
// surface kinetics mechanism
|
||||
//
|
||||
// HKM Should do it here for all kinetics objects so that
|
||||
// bulk will eventually work.
|
||||
if (do_time) {
|
||||
kindexSP = 0;
|
||||
for (isp = 0; isp < m_numSurfPhases; isp++) {
|
||||
|
|
@ -546,13 +496,9 @@ void solveSP::resjac_eval(SquareMatrix& jac,
|
|||
{
|
||||
size_t kColIndex = 0, nsp, jsp, i, kCol;
|
||||
doublereal dc, cSave, sd;
|
||||
/*
|
||||
* Calculate the residual
|
||||
*/
|
||||
// Calculate the residual
|
||||
fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
|
||||
/*
|
||||
* Now we will look over the columns perturbing each unknown.
|
||||
*/
|
||||
// Now we will look over the columns perturbing each unknown.
|
||||
for (jsp = 0; jsp < m_numSurfPhases; jsp++) {
|
||||
nsp = m_nSpeciesSurfPhase[jsp];
|
||||
sd = m_ptrsSurfPhase[jsp]->siteDensity();
|
||||
|
|
@ -607,16 +553,12 @@ static doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, i
|
|||
*label = -1;
|
||||
|
||||
for (size_t i = 0; i < dim; i++) {
|
||||
/*
|
||||
* Calculate the new suggested new value of x[i]
|
||||
*/
|
||||
// Calculate the new suggested new value of x[i]
|
||||
xnew = x[i] - damp * dxneg[i];
|
||||
|
||||
/*
|
||||
* Calculate the allowed maximum and minimum values of x[i]
|
||||
* - Only going to allow x[i] to converge to zero by a
|
||||
* single order of magnitude at a time
|
||||
*/
|
||||
// Calculate the allowed maximum and minimum values of x[i]
|
||||
// - Only going to allow x[i] to converge to zero by a
|
||||
// single order of magnitude at a time
|
||||
xtop = 1.0 - 0.1*fabs(1.0-x[i]);
|
||||
xbot = fabs(x[i]*0.1) - 1.0e-16;
|
||||
if (xnew > xtop) {
|
||||
|
|
@ -631,27 +573,23 @@ static doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, i
|
|||
}
|
||||
}
|
||||
damp = std::max(damp, 1e-2);
|
||||
/*
|
||||
* Only allow the damping parameter to increase by a factor of three each
|
||||
* iteration. Heuristic to avoid oscillations in the value of damp
|
||||
*/
|
||||
|
||||
// Only allow the damping parameter to increase by a factor of three each
|
||||
// iteration. Heuristic to avoid oscillations in the value of damp
|
||||
if (damp > damp_old*3) {
|
||||
damp = damp_old*3;
|
||||
*label = -1;
|
||||
}
|
||||
|
||||
/*
|
||||
* Save old value of the damping parameter for use
|
||||
* in subsequent calls.
|
||||
*/
|
||||
// Save old value of the damping parameter for use in subsequent calls.
|
||||
damp_old = damp;
|
||||
return damp;
|
||||
|
||||
} /* calc_damping */
|
||||
|
||||
/*
|
||||
* This function calculates the norm of an update, dx[],
|
||||
* based on the weighted values of x.
|
||||
* This function calculates the norm of an update, dx[], based on the
|
||||
* weighted values of x.
|
||||
*/
|
||||
static doublereal calcWeightedNorm(const doublereal wtX[], const doublereal dx[], size_t dim)
|
||||
{
|
||||
|
|
@ -673,10 +611,9 @@ void solveSP::calcWeights(doublereal wtSpecies[], doublereal wtResid[],
|
|||
{
|
||||
size_t k, jcol, kindex, isp, nsp;
|
||||
doublereal sd;
|
||||
/*
|
||||
* First calculate the weighting factor for the concentrations of
|
||||
* the surface species and bulk species.
|
||||
*/
|
||||
|
||||
// First calculate the weighting factor for the concentrations of the
|
||||
// surface species and bulk species.
|
||||
kindex = 0;
|
||||
for (isp = 0; isp < m_numSurfPhases; isp++) {
|
||||
nsp = m_nSpeciesSurfPhase[isp];
|
||||
|
|
@ -694,12 +631,10 @@ void solveSP::calcWeights(doublereal wtSpecies[], doublereal wtResid[],
|
|||
}
|
||||
}
|
||||
}
|
||||
/*
|
||||
* Now do the residual Weights. Since we have the Jacobian, we
|
||||
* will use it to generate a number based on the what a significant
|
||||
* change in a solution variable does to each residual.
|
||||
* This is a row sum scale operation.
|
||||
*/
|
||||
|
||||
// Now do the residual Weights. Since we have the Jacobian, we will use it
|
||||
// to generate a number based on the what a significant change in a solution
|
||||
// variable does to each residual. This is a row sum scale operation.
|
||||
for (k = 0; k < m_neq; k++) {
|
||||
wtResid[k] = 0.0;
|
||||
for (jcol = 0; jcol < m_neq; jcol++) {
|
||||
|
|
@ -752,10 +687,8 @@ doublereal solveSP::calc_t(doublereal netProdRateSolnSP[],
|
|||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* Increase time step exponentially as same species repeatedly
|
||||
* controls time step
|
||||
*/
|
||||
// Increase time step exponentially as same species repeatedly controls time
|
||||
// step
|
||||
if (*label == *label_old) {
|
||||
*label_factor *= 1.5;
|
||||
} else {
|
||||
|
|
@ -763,7 +696,7 @@ doublereal solveSP::calc_t(doublereal netProdRateSolnSP[],
|
|||
*label_factor = 1.0;
|
||||
}
|
||||
return inv_timeScale / *label_factor;
|
||||
} /* calc_t */
|
||||
} // calc_t
|
||||
|
||||
void solveSP::print_header(int ioflag, int ifunc, doublereal time_scale,
|
||||
int damping, doublereal reltol, doublereal abstol)
|
||||
|
|
@ -859,6 +792,6 @@ void solveSP::printIteration(int ioflag, doublereal damp, int label_d,
|
|||
}
|
||||
writelog("\n");
|
||||
}
|
||||
} /* printIteration */
|
||||
} // printIteration
|
||||
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue