From 651ce785cc6ccf8a3975f0ff28f8a6f0f1c279f7 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 18 Apr 2013 22:07:43 +0000 Subject: [PATCH] Cleaned up Doxygen docs for class ChemEquil --- include/cantera/equil/ChemEquil.h | 126 ++++++++++++++++++++++++----- src/equil/ChemEquil.cpp | 129 ------------------------------ 2 files changed, 106 insertions(+), 149 deletions(-) diff --git a/include/cantera/equil/ChemEquil.h b/include/cantera/equil/ChemEquil.h index 1d10f533e..3e434996a 100644 --- a/include/cantera/equil/ChemEquil.h +++ b/include/cantera/equil/ChemEquil.h @@ -1,17 +1,13 @@ /** - * @file ChemEquil.h - * - * Chemical equilibrium. + * @file ChemEquil.h Chemical equilibrium. */ /* * Copyright 2001 California Institute of Technology */ - #ifndef CT_CHEM_EQUIL_H #define CT_CHEM_EQUIL_H - // Cantera includes #include "cantera/base/ct_defs.h" #include "cantera/base/vec_functions.h" @@ -26,6 +22,7 @@ namespace Cantera { +/// map property strings to integers int _equilflag(const char* xy); /** @@ -97,9 +94,7 @@ class PropertyCalculator; */ class ChemEquil { - public: - //! Default Constructor ChemEquil(); //! Constructor combined with the initialization function @@ -112,8 +107,31 @@ public: virtual ~ChemEquil(); + /*! + * Equilibrate a phase, holding the elemental composition fixed + * at the initial value found within the ThermoPhase object *s*. + * + * The value of 2 specified properties are obtained by querying the + * ThermoPhase object. The properties must be already contained + * within the current thermodynamic state of the system. + */ int equilibrate(thermo_t& s, const char* XY, bool useThermoPhaseElementPotentials = false, int loglevel = 0); + + /*! + * Compute the equilibrium composition for 2 specified + * properties and the specified element moles. + * + * The 2 specified properties are obtained by querying the + * ThermoPhase object. The properties must be already contained + * within the current thermodynamic state of the system. + * + * @param elMoles specified vector of element abundances. + * + * @return Successful returns are indicated by a return value of 0. + * Unsuccessful returns are indicated by a return value of -1 for lack + * of convergence or -3 for a singular jacobian. + */ int equilibrate(thermo_t& s, const char* XY, vector_fp& elMoles, bool useThermoPhaseElementPotentials = false, int loglevel = 0); const vector_fp& elementPotentials() const { @@ -139,27 +157,91 @@ protected: */ thermo_t* m_phase; - /// number of atoms of element m in species k. + //! number of atoms of element m in species k. doublereal nAtoms(size_t k, size_t m) const { return m_comp[k*m_mm + m]; } + /*! + * Prepare for equilibrium calculations. + * @param s object representing the solution phase. + */ void initialize(thermo_t& s); + /*! + * Set mixture to an equilibrium state consistent with specified + * element potentials and temperature. + * + * @param x vector of non-dimensional element potentials + * \f[ \lambda_m/RT \f]. + * @param t temperature in K. + */ void setToEquilState(thermo_t& s, const vector_fp& x, doublereal t); + //! Estimate the initial mole numbers. This version borrows from the + //! MultiPhaseEquil solver. int setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, int loglevel = 0); + //! Generate a starting estimate for the element potentials. int estimateElementPotentials(thermo_t& s, vector_fp& lambda, vector_fp& elMolesGoal, int loglevel = 0); + /*! + * Do a calculation of the element potentials using the Brinkley method, + * p. 129 Smith and Missen. + * + * We have found that the previous estimate may not be good enough to + * avoid drastic numerical issues associated with the use of a numerically + * generated jacobian used in the main algorithm. + * + * The Brinkley algorithm, here, assumes a constant T, P system and uses a + * linearized analytical Jacobian that turns out to be very stable even + * given bad initial guesses. + * + * The pressure and temperature to be used are in the ThermoPhase object + * input into the routine. + * + * The initial guess for the element potentials used by this routine is + * taken from the input vector, x. + * + * elMoles is the input element abundance vector to be matched. + * + * Nonideal phases are handled in principle. This is done by calculating + * the activity coefficients and adding them into the formula in the + * correct position. However, these are treated as a rhs contribution + * only. Therefore, convergence might be a problem. This has not been + * tested. Also molality based unit systems aren't handled. + * + * On return, int return value contains the success code: + * - 0 - successful + * - 1 - unsuccessful, max num iterations exceeded + * - -3 - unsuccessful, singular jacobian + * + * NOTE: update for activity coefficients. + */ int estimateEP_Brinkley(thermo_t& s, vector_fp& lambda, vector_fp& elMoles); + //! Find an acceptable step size and take it. + /*! + * The original implementation employed a line search technique that + * enforced a reduction in the norm of the residual at every successful + * step. Unfortunately, this method created false convergence errors near + * the end of a significant number of steps, usually special conditions + * where there were stoichiometric constraints. + * + * This new method just does a delta damping approach, based on limiting + * the jump in the dimensionless element potentials. Mole fractions are + * limited to a factor of 2 jump in the values from this method. Near + * convergence, the delta damping gets out of the way. + */ int dampStep(thermo_t& s, vector_fp& oldx, double oldf, vector_fp& grad, vector_fp& step, vector_fp& x, double& f, vector_fp& elmols, double xval, double yval); + /** + * Evaluates the residual vector F, of length #m_mm + */ void equilResidual(thermo_t& s, const vector_fp& x, const vector_fp& elmtotal, vector_fp& resid, double xval, double yval, int loglevel = 0); @@ -170,20 +252,27 @@ protected: void adjustEloc(thermo_t& s, vector_fp& elMolesGoal); + //! Update internally stored state information. void update(const thermo_t& s); + /** + * Given a vector of dimensionless element abundances, this routine + * calculates the moles of the elements and the moles of the species. + * + * @param[in] x = current dimensionless element potentials.. + */ double calcEmoles(thermo_t& s, vector_fp& x, const double& n_t, const vector_fp& Xmol_i_calc, vector_fp& eMolesCalc, vector_fp& n_i_calc, double pressureConst); - size_t m_mm; - size_t m_kk; + size_t m_mm; //!< number of elements in the phase + size_t m_kk; //!< number of species in the phase size_t m_skip; /** * This is equal to the rank of the stoichiometric coefficient - * matrix when it is computed. It's initialized to m_mm. + * matrix when it is computed. It's initialized to #m_mm. */ size_t m_nComponents; @@ -191,12 +280,12 @@ protected: /** * Current value of the mole fractions in the single phase. - * -> length = m_kk. + * -> length = #m_kk. */ vector_fp m_molefractions; /** * Current value of the dimensional element potentials - * -> length = m_mm + * -> length = #m_mm */ vector_fp m_lambda; @@ -213,6 +302,7 @@ protected: vector_fp m_reswork; vector_fp m_jwork1; vector_fp m_jwork2; + /* * Storage of the element compositions * natom(k,m) = m_comp[k*m_mm+ m]; @@ -220,6 +310,7 @@ protected: vector_fp m_comp; doublereal m_temp, m_dens; doublereal m_p0; + /** * Index of the element id corresponding to the electric charge of each * species. Equal to -1 if there is no such element id. @@ -230,6 +321,7 @@ protected: vector_fp m_grt; vector_fp m_mu_RT; + /** * Dimensionless values of the gibbs free energy for the * standard state of each species, at the temperature and @@ -238,18 +330,12 @@ protected: vector_fp m_muSS_RT; std::vector m_component; - /* - * element fractional cutoff, below which the element will be - * zeroed. - */ + //! element fractional cutoff, below which the element will be zeroed. double m_elemFracCutoff; bool m_doResPerturb; - std::vector m_orderVectorElements; std::vector m_orderVectorSpecies; - - }; extern int ChemEquil_print_lvl; diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index 627eddbe6..e2dca056d 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -33,7 +33,6 @@ int Cantera::ChemEquil_print_lvl = 0; namespace Cantera { -/// map property strings to integers int _equilflag(const char* xy) { string flag = string(xy); @@ -57,25 +56,12 @@ int _equilflag(const char* xy) return -1; } - -//----------------------------------------------------------- -// construction / destruction -//----------------------------------------------------------- - - -/// Default Constructor. ChemEquil::ChemEquil() : m_skip(npos), m_elementTotalSum(1.0), m_p0(OneAtm), m_eloc(npos), m_elemFracCutoff(1.0E-100), m_doResPerturb(false) {} -//! Constructor combined with the initialization function -/*! - * This constructor initializes the ChemEquil object with everything it - * needs to start solving equilibrium problems. - * @param s ThermoPhase object that will be used in the equilibrium calls. - */ ChemEquil::ChemEquil(thermo_t& s) : m_skip(npos), m_elementTotalSum(1.0), @@ -90,10 +76,6 @@ ChemEquil::~ChemEquil() { } -/** - * Prepare for equilibrium calculations. - * @param s object representing the solution phase. - */ void ChemEquil::initialize(thermo_t& s) { // store a pointer to s and some of its properties locally. @@ -169,16 +151,6 @@ void ChemEquil::initialize(thermo_t& s) } } - -/** - * Set mixture to an equilibrium state consistent with specified - * element potentials and temperature. - * - * @param lambda_RT vector of non-dimensional element potentials - * \f[ \lambda_m/RT \f]. - * @param t temperature in K. - * - */ void ChemEquil::setToEquilState(thermo_t& s, const vector_fp& lambda_RT, doublereal t) { @@ -199,10 +171,6 @@ void ChemEquil::setToEquilState(thermo_t& s, update(s); } - -/** - * update internally stored state information. - */ void ChemEquil::update(const thermo_t& s) { @@ -233,8 +201,6 @@ void ChemEquil::update(const thermo_t& s) } } -/// Estimate the initial mole numbers. This version borrows from the -/// MultiPhaseEquil solver. int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, int loglevel) { @@ -302,10 +268,6 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, return iok; } - -/** - * Generate a starting estimate for the element potentials. - */ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, vector_fp& elMolesGoal, int loglevel) { @@ -440,15 +402,6 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, return info; } - -/** - * Equilibrate a phase, holding the elemental composition fixed - * at the initial value found within the ThermoPhase object. - * - * The value of 2 specified properties are obtained by querying the - * ThermoPhase object. The properties must be already contained - * within the current thermodynamic state of the system. - */ int ChemEquil::equilibrate(thermo_t& s, const char* XY, bool useThermoPhaseElementPotentials, int loglevel) { @@ -461,22 +414,6 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XY, loglevel-1); } - -/** - * Compute the equilibrium composition for 2 specified - * properties and the specified element moles. - * - * elMoles = specified vector of element abundances. - * - * The 2 specified properties are obtained by querying the - * ThermoPhase object. The properties must be already contained - * within the current thermodynamic state of the system. - * - * Return variable: - * Successful returns are indicated by a return value of 0. - * Unsuccessful returns are indicated by a return value of -1 for - * lack of convergence or -3 for a singular jacobian. - */ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, vector_fp& elMolesGoal, bool useThermoPhaseElementPotentials, @@ -1053,20 +990,6 @@ converge: goto next; } - -/* - * dampStep: Come up with an acceptable step size. The original implementation - * employed a line search technique that enforced a reduction in the - * norm of the residual at every successful step. Unfortunately, - * this method created false convergence errors near the end of - * a significant number of steps, usually special conditions where - * there were stoichiometric constraints. - * - * This new method just does a delta damping approach, based on limiting - * the jump in the dimensionless element potentials. Mole fractions are - * limited to a factor of 2 jump in the values from this method. - * Near convergence, the delta damping gets out of the way. - */ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx, double oldf, vector_fp& grad, vector_fp& step, vector_fp& x, double& f, vector_fp& elmols, double xval, double yval) @@ -1113,10 +1036,6 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx, return 1; } - -/** - * Evaluates the residual vector F, of length mm - */ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, const vector_fp& elmFracGoal, vector_fp& resid, doublereal xval, doublereal yval, int loglevel) @@ -1189,9 +1108,6 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, #endif } - -//-------------------- Jacobian evaluation --------------------------- - void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, const vector_fp& elmols, DenseMatrix& jac, doublereal xval, doublereal yval, int loglevel) @@ -1239,14 +1155,6 @@ void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, } } -/** - * Given a vector of dimensionless element abundances, - * this routine calculates the moles of the elements and - * the moles of the species. - * Input - * -------- - * x[m] = current dimensionless element potentials.. - */ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t, const vector_fp& Xmol_i_calc, vector_fp& eMolesCalc, vector_fp& n_i_calc, @@ -1287,42 +1195,6 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t, return n_t_calc; } -/** - * Do a calculation of the element potentials using - * the Brinkley method, p. 129 Smith and Missen. - * - * We have found that the previous estimate may not be good - * enough to avoid drastic numerical issues associated with - * the use of a numerically generated jacobian used in the - * main algorithm. - * - * The Brinkley algorithm, here, assumes a constant T, P system - * and uses a linearized analytical Jacobian that turns out - * to be very stable even given bad initial guesses. - * - * The pressure and temperature to be used are in the - * ThermoPhase object input into the routine. - * - * The initial guess for the element potentials - * used by this routine is taken from the - * input vector, x. - * - * elMoles is the input element abundance vector to be matched. - * - * Nonideal phases are handled in principle. This is done by - * calculating the activity coefficients and adding them - * into the formula in the correct position. However, - * these are treated as a rhs contribution only. Therefore, - * convergence might be a problem. This has not been tested. - * Also molality based unit systems aren't handled. - * - * On return, int return value contains the success code: - * 0 - successful - * 1 - unsuccessful, max num iterations exceeded - * -3 - unsuccessful, singular jacobian - * - * NOTE: update for activity coefficients. - */ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, vector_fp& elMoles) { @@ -1998,4 +1870,3 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal) } } // namespace -