diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index daf20d353..45c1801e3 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -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 diff --git a/include/cantera/kinetics/ImplicitSurfChem.h b/include/cantera/kinetics/ImplicitSurfChem.h index 87a28ca52..f2cc1757d 100644 --- a/include/cantera/kinetics/ImplicitSurfChem.h +++ b/include/cantera/kinetics/ImplicitSurfChem.h @@ -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 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 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: diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index b3d65b7cd..89b9537db 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -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. diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index a1aa75877..c7e0dee70 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -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 "" 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 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 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 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; diff --git a/include/cantera/kinetics/RateCoeffMgr.h b/include/cantera/kinetics/RateCoeffMgr.h index da3b97098..16988de27 100644 --- a/include/cantera/kinetics/RateCoeffMgr.h +++ b/include/cantera/kinetics/RateCoeffMgr.h @@ -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; diff --git a/include/cantera/kinetics/ReactionPath.h b/include/cantera/kinetics/ReactionPath.h index 94e08e4a2..835b14c8e 100644 --- a/include/cantera/kinetics/ReactionPath.h +++ b/include/cantera/kinetics/ReactionPath.h @@ -85,8 +85,7 @@ public: typedef std::map 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); diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index 413985513..358feb94a 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -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); diff --git a/include/cantera/kinetics/StoichManager.h b/include/cantera/kinetics/StoichManager.h index b70284db1..855d81599 100644 --- a/include/cantera/kinetics/StoichManager.h +++ b/include/cantera/kinetics/StoichManager.h @@ -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 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& 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& k, const vector_fp& order, const vector_fp& stoich) { diff --git a/include/cantera/kinetics/importKinetics.h b/include/cantera/kinetics/importKinetics.h index ae2155896..229c871d8 100644 --- a/include/cantera/kinetics/importKinetics.h +++ b/include/cantera/kinetics/importKinetics.h @@ -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 diff --git a/include/cantera/kinetics/reaction_defs.h b/include/cantera/kinetics/reaction_defs.h index f90dacf93..5767d27b7 100644 --- a/include/cantera/kinetics/reaction_defs.h +++ b/include/cantera/kinetics/reaction_defs.h @@ -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; diff --git a/include/cantera/kinetics/solveSP.h b/include/cantera/kinetics/solveSP.h index 3b06f4b83..692141b39 100644 --- a/include/cantera/kinetics/solveSP.h +++ b/include/cantera/kinetics/solveSP.h @@ -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 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) */ diff --git a/src/kinetics/AqueousKinetics.cpp b/src/kinetics/AqueousKinetics.cpp index f3617d96b..1cb4a74a0 100644 --- a/src/kinetics/AqueousKinetics.cpp +++ b/src/kinetics/AqueousKinetics.cpp @@ -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 diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index 2f947f27b..32123ee41 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -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) { diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 46dbc62f8..0354dab4c 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -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 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); diff --git a/src/kinetics/Group.cpp b/src/kinetics/Group.cpp index 0d51f01a4..a5f8bcf86 100644 --- a/src/kinetics/Group.cpp +++ b/src/kinetics/Group.cpp @@ -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 diff --git a/src/kinetics/ImplicitSurfChem.cpp b/src/kinetics/ImplicitSurfChem.cpp index 3d94190ec..897124f61 100644 --- a/src/kinetics/ImplicitSurfChem.cpp +++ b/src/kinetics/ImplicitSurfChem.cpp @@ -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) { diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 5e6f7358b..6e84bc494 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -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(nReactions(), 1); diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index b545d0ef2..4536b9d46 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -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; } diff --git a/src/kinetics/KineticsFactory.cpp b/src/kinetics/KineticsFactory.cpp index 3a9534f0e..ebcc44688 100644 --- a/src/kinetics/KineticsFactory.cpp +++ b/src/kinetics/KineticsFactory.cpp @@ -22,17 +22,15 @@ std::mutex KineticsFactory::kinetics_mutex; Kinetics* KineticsFactory::newKinetics(XML_Node& phaseData, vector 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 diff --git a/src/kinetics/ReactionPath.cpp b/src/kinetics/ReactionPath.cpp index e8602fed7..e6717db61 100644 --- a/src/kinetics/ReactionPath.cpp +++ b/src/kinetics/ReactionPath.cpp @@ -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 > allProducts(m_nr); vector > 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 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 >& 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); } diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index d25bba29e..58f6b398d 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -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 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 incl = rxns.getChildren("include"); vector 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 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 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 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 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; } diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp index 55b120958..fd5cda06d 100644 --- a/src/kinetics/solveSP.cpp +++ b/src/kinetics/solveSP.cpp @@ -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 }