diff --git a/Cantera/src/IdealGasPhase.h b/Cantera/src/IdealGasPhase.h index 0b63a0d51..018ed2bd6 100644 --- a/Cantera/src/IdealGasPhase.h +++ b/Cantera/src/IdealGasPhase.h @@ -1,7 +1,5 @@ /** - * - * @file IdealGasPhase.h - * ` + * @file IdealGasPhase.h * ThermoPhase object for the ideal gas equation of state. */ @@ -26,15 +24,18 @@ namespace Cantera { - //!Class %IdealGasPhase represents low-density gases that obey the + //! Class %IdealGasPhase represents low-density gases that obey the //! ideal gas equation of state. /*! * * %IdealGasPhase derives from class ThermoPhase, * and overloads the virtual methods defined there with ones that * use expressions appropriate for ideal gas mixtures. - * - * This class is optimized for speed of execution. + * + * The independent unknowns are density, mass fraction, and temperature. + * the #setPressure() function will calculate the density consistent with + * the current mass fraction vector and temperature and the desired pressure, + * and then set the density in the derived State object. * *
*

Specification of Species Standard %State Properties

@@ -47,9 +48,16 @@ namespace Cantera { * description of the specification of reference state species thermodynamics functions). * The reference state, * where the pressure is fixed at a single pressure, - * is key species property calculation for the Ideal Gas Equation + * is a key species property calculation for the Ideal Gas Equation * of state. * + * This class is optimized for speed of execution. All calls to thermodynamic functions + * first call internal routines (aka #enthalpy_RT_ref()) which return references + * the reference state thermodynamics functions. Within these internal reference + * state functions, the function #_updateThermo() is called, that first checks to see + * whether the temperature has changed. If it has, it updates the internal reference + * state thermo functions by calling the SpeciesThermo object. + * * Functions for the calculation of standard state properties for species * at arbitray pressure are provided in %IdealGasPhase. However, they * are all derived from their reference state conterparts. @@ -87,10 +95,10 @@ namespace Cantera { * \mu^{ref}_k(T) = h^{ref}_k(T) - T S^{ref}_k(T) * \f] * - * The standard state internal energy is obtained from the enthalpy function too + * The standard state internal energy is obtained from the enthalpy function also * * \f[ - * u^o_k(T,P) = h^o_k(T) - R T + * u^o_k(T,P) = h^o_k(T) - R T * \f] * * The molar volume of a species is given by the ideal gas law @@ -124,19 +132,25 @@ namespace Cantera { * \mu_k(T,P) = \mu^{ref}_k(T, P) + R T \log(\frac{P X_k}{P_{ref}}) * \f] * - * The partial molar entropy for species k is given by the following relation, + * The partial molar entropy for species k is given by the following relation, * * \f[ * \tilde{s}_k(T,P) = s^o_k(T,P) - R \log(X_k) = s^{ref}_k(T) - R \log(\frac{P X_k}{P_{ref}}) * \f] * - * The partial molar enthalpy for species k is + * The partial molar enthalpy for species k is * * \f[ * \tilde{h}_k(T,P) = h^o_k(T,P) = h^{ref}_k(T) * \f] * - * The partial molar heat capacity for species k is + * The partial molar Internal Energy for species k is + * + * \f[ + * \tilde{u}_k(T,P) = u^o_k(T,P) = u^{ref}_k(T) + * \f] + * + * The partial molar Heat Capacity for species k is * * \f[ * \tilde{Cp}_k(T,P) = Cp^o_k(T,P) = Cp^{ref}_k(T) @@ -146,7 +160,7 @@ namespace Cantera { *
*

%Application within %Kinetics Managers

*
- + * * \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k / * C^s_k, \f$ where \f$ C^s_k \f$ is a standard concentration * defined below and \f$ a_k \f$ are activities used in the @@ -246,12 +260,49 @@ namespace Cantera { *

Instantiation of the Class

*
* + * + * The constructor for this phase is located in the default ThermoFactory + * for %Cantera. A new %IdealGasPhase may be created by the following code snippet: * + * @code + * XML_Node * const xs = xc->findNameID("phase", "silane"); + * ThermoPhase *silane_tp = newPhase(*xs); + * IdealGasPhase *silaneGas = dynamic_cast (silane_tp); + * @endcode + * + * or by the following constructor: + * + * @code + * XML_Node * const xs = xc->findNameID("phase", "silane"); + * IdealGasPhase *silaneGas = new IdealGasPhase(*xs); + * @endcode + *
*

XML Example

*
+ * An example of an XML Element named phase setting up a IdealGasPhase object named silane + * is given below. + * + * @verbatim + + + Si H He + + H2 H HE SIH4 SI SIH SIH2 SIH3 H3SISIH SI2H6 + H2SISIH2 SI3H8 SI2 SI3 + + + + + + + @endverbatim + * + * The model attribute "IdealGas" of the thermo XML element identifies the phase as + * being of the type handled by the IdealGasPhase object. + * + * @ingroup thermoprops * - * @ingroup thermoprops */ class IdealGasPhase : public ThermoPhase { @@ -698,19 +749,34 @@ namespace Cantera { virtual void getStandardVolumes_ref(doublereal *vol) const; //@} - /// @name New Methods Defined Here ------------------------------------------------- + /// @name NonVirtual Internal methods to Return References to Reference State Thermo //@{ + //! Returns a reference to the dimensionless reference state enthalpy vector. + /*! + * This function is part of the layer that checks/recalculates the reference + * state thermo functions. + */ const array_fp& enthalpy_RT_ref() const { _updateThermo(); return m_h0_RT; } + //! Returns a reference to the dimensionless reference state Gibbs free energy vector. + /*! + * This function is part of the layer that checks/recalculates the reference + * state thermo functions. + */ const array_fp& gibbs_RT_ref() const { _updateThermo(); return m_g0_RT; } + //! Returns a reference to the exponent of the dimensionless reference state Gibbs Free energy vector. + /*! + * This function is part of the layer that checks/recalculates the reference + * state thermo functions. + */ const array_fp& expGibbs_RT_ref() const { _updateThermo(); int k; @@ -718,26 +784,42 @@ namespace Cantera { return m_expg0_RT; } + //! Returns a reference to the dimensionless reference state Entropy vector. + /*! + * This function is part of the layer that checks/recalculates the reference + * state thermo functions. + */ const array_fp& entropy_R_ref() const { _updateThermo(); return m_s0_R; } + //! Returns a reference to the dimensionless reference state Heat Capacity vector. + /*! + * This function is part of the layer that checks/recalculates the reference + * state thermo functions. + */ const array_fp& cp_R_ref() const { _updateThermo(); return m_cp0_R; } - // @} + //@} - /** - * @internal Initialize. This method is provided to allow + //! Initialize the ThermoPhase object after all species have been set up + /*! + * @internal Initialize. + * + * This method is provided to allow * subclasses to perform any initialization required after all * species have been added. For example, it might be used to * resize internal work arrays that must have an entry for * each species. The base class implementation does nothing, * and subclasses that do not require initialization do not - * need to overload this method. + * need to overload this method. When importing a CTML phase + * description, this method is called from ThermoPhase::initThermoXML(), + * which is called from importPhase(), + * just prior to returning from function importPhase(). * * @see importCTML.cpp */ @@ -829,7 +911,14 @@ namespace Cantera { private: + //! Update the species reference state thermodynamic functions + /*! + * The polynomials for the standard state functions are only + * reevalulated if the temperature has changed. + * + */ void _updateThermo() const; + }; } diff --git a/Cantera/src/SurfPhase.h b/Cantera/src/SurfPhase.h index a91993b88..bf5c49b64 100644 --- a/Cantera/src/SurfPhase.h +++ b/Cantera/src/SurfPhase.h @@ -537,7 +537,7 @@ namespace Cantera { private: - //! Update the species standard state thermodynamic functions + //! Update the species reference state thermodynamic functions /*! * The polynomials for the standard state functions are only * reevalulated if the temperature has changed. diff --git a/Cantera/src/equil.h b/Cantera/src/equil.h index 45ee613af..cd43ba7c5 100644 --- a/Cantera/src/equil.h +++ b/Cantera/src/equil.h @@ -6,7 +6,7 @@ ***********************************************************************/ // Copyright 2001 California Institute of Technology - /** +/** * @file equil.h * This file contains the definition of some high level general equilibration * routines and the text for the module \ref equilfunctions. @@ -22,27 +22,84 @@ namespace Cantera { - /*! - * @defgroup equilfunctions Equilibrium Solver Capability - * - * Cantera has several different equilibrium routines. - */ - //----------------------------------------------------------- - // convenience functions - //----------------------------------------------------------- - - /** - * Return variable is equal to the number of subroutine attempts - * it took to equilibrate the system. - */ - int equilibrate(thermo_t& s, const char* XY, - int solver = -1, doublereal rtol = 1.0e-9, int maxsteps = 1000, - int maxiter = 100, int loglevel = -99); - - doublereal equilibrate(MultiPhase& s, const char* XY, - doublereal tol = 1.0e-9, int maxsteps = 1000, int maxiter = 100, - int loglevel = -99); + /*! + * @defgroup equilfunctions Equilibrium Solver Capability + * + * Cantera has several different equilibrium routines. + */ + //----------------------------------------------------------- + // convenience functions + //----------------------------------------------------------- + + //! Equilibrate a ThermoPhase object + /*! + * Set a single-phase chemical solution to chemical equilibrium. + * This is a convenience function that uses one or the other of + * the two chemical equilibrium solvers. The XY parameter indicates what two + * thermodynamic quantities, other than element composition, are to be held + * constant during the equilibration process. + * + * @param s ThermoPhase object that will be equilibrated. + * @param XY String representation of what two properties + * are being held constant + * @param solver ID of the solver to be used to equlibrate the phase. + * If solver = 0, the ChemEquil solver will be used, + * and if solver = 1, the + * MultiPhaseEquil solver will be used (slower than ChemEquil, + * but more stable). If solver < 0 (default, then ChemEquil will + * be tried first, and if it fails MultiPhaseEquil will be tried. + * @param rtol Relative tolerance + * @param maxsteps Maximum number of steps to take to find the solution + * @param maxiter For the MultiPhaseEquil solver only, this is + * the maximum number of outer temperature or pressure iterations + * to take when T and/or P is not held fixed. + * @param loglevel loglevel Controls amount of diagnostic output. loglevel + * = 0 suppresses diagnostics, and increasingly-verbose messages + * are written as loglevel increases. The messages are written to + * a file in HTML format for viewing in a web browser. + * @see HTML_logs + * + * @return + * Return variable is equal to the number of subroutine attempts + * it took to equilibrate the system. + * + * + * @ingroup equilfunctions + * @ingroup equil + */ + int equilibrate(thermo_t& s, const char* XY, + int solver = -1, doublereal rtol = 1.0e-9, int maxsteps = 1000, + int maxiter = 100, int loglevel = -99); + //! Equilibrate a MultiPhase object + /*! + * Equilibrate a MultiPhase object. The XY parameter indicates what two + * thermodynamic quantities, other than element composition, are to be held + * constant during the equilibration process. + * + * This is the top-level driver for multiphase equilibrium. It + * doesn't do much more than call the equilibrate method of class + * MultiPhase, except that it adds some messages to the logfile, + * if loglevel is set > 0. + * + * @param s MultiPhase object that will be equilibrated. + * @param XY String representation of what is being held constant + * @param rtol Relative tolerance + * @param maxsteps Maximum number of steps + * @param maxiter Maximum iterations + * @param loglevel loglevel + * + * @return + * Return variable is equal to the number of subroutine attempts + * it took to equilibrate the system. + * + * @ingroup equilfunctions + * @ingroup equil + */ + doublereal equilibrate(MultiPhase& s, const char* XY, + doublereal rtol = 1.0e-9, int maxsteps = 1000, int maxiter = 100, + int loglevel = -99); + } #endif diff --git a/Cantera/src/equilibrate.cpp b/Cantera/src/equilibrate.cpp index 2e383d3d1..e0769d809 100644 --- a/Cantera/src/equilibrate.cpp +++ b/Cantera/src/equilibrate.cpp @@ -1,6 +1,5 @@ /** * @file equilibrate.cpp - * * Driver routines for the chemical equilibrium solvers. * */ @@ -12,177 +11,177 @@ namespace Cantera { - /** - * Set a multiphase mixture to a state of chemical equilibrium. - * This is the top-level driver for multiphase equilibrium. It - * doesn't do much more than call the equilibrate method of class - * MultiPhase, except that it adds some messages to the logfile, - * if loglevel is set > 0. - * - * @ingroup equil - */ - doublereal equilibrate(MultiPhase& s, const char* XY, - doublereal tol, int maxsteps, int maxiter, - int loglevel) { + /* + * Set a multiphase mixture to a state of chemical equilibrium. + * This is the top-level driver for multiphase equilibrium. It + * doesn't do much more than call the equilibrate method of class + * MultiPhase, except that it adds some messages to the logfile, + * if loglevel is set > 0. + * + * @ingroup equil + */ + doublereal equilibrate(MultiPhase& s, const char* XY, + doublereal tol, int maxsteps, int maxiter, + int loglevel) { - beginLogGroup("equilibrate",loglevel); - addLogEntry("multiphase equilibrate function"); - beginLogGroup("arguments"); - addLogEntry("XY",XY); - addLogEntry("tol",tol); - addLogEntry("maxsteps",maxsteps); - addLogEntry("maxiter",maxiter); - addLogEntry("loglevel",loglevel); - endLogGroup("arguments"); + beginLogGroup("equilibrate",loglevel); + addLogEntry("multiphase equilibrate function"); + beginLogGroup("arguments"); + addLogEntry("XY",XY); + addLogEntry("tol",tol); + addLogEntry("maxsteps",maxsteps); + addLogEntry("maxiter",maxiter); + addLogEntry("loglevel",loglevel); + endLogGroup("arguments"); - s.init(); - int ixy = _equilflag(XY); - if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) { - try { - double err = s.equilibrate(ixy, tol, maxsteps, maxiter); - addLogEntry("Success. Error",err); - endLogGroup("equilibrate"); - return err; - } - catch (CanteraError e) { - addLogEntry("Failure.",lastErrorMessage()); - endLogGroup("equilibrate"); - throw e; - } - } - else { - addLogEntry("multiphase equilibrium can be done only for TP, HP, SP, or TV"); - endLogGroup("equilibrate"); - throw CanteraError("equilibrate","unsupported option"); - return -1.0; - } + s.init(); + int ixy = _equilflag(XY); + if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) { + try { + double err = s.equilibrate(ixy, tol, maxsteps, maxiter); + addLogEntry("Success. Error",err); + endLogGroup("equilibrate"); + return err; + } + catch (CanteraError e) { + addLogEntry("Failure.",lastErrorMessage()); + endLogGroup("equilibrate"); + throw e; + } } - - - /// Set a single-phase chemical solution to chemical equilibrium. - /// This is a convenience function that uses one or the other of - /// the two chemical equilibrium solvers. - /// - /// @param s The object to set to an equilibrium state - /// - /// @param XY An integer specifying the two properties to be held - /// constant. - /// - /// @param solver The equilibrium solver to use. If solver = 0, - /// the ChemEquil solver will be used, and if solver = 1, the - /// MultiPhaseEquil solver will be used (slower than ChemEquil, - /// but more stable). If solver < 0 (default, then ChemEquil will - /// be tried first, and if it fails MultiPhaseEquil will be tried. - /// - /// @param maxsteps The maximum number of steps to take to find - /// the solution. - /// - /// @param maxiter For the MultiPhaseEquil solver only, this is - /// the maximum number of outer temperature or pressure iterations - /// to take when T and/or P is not held fixed. - /// - /// @param loglevel Controls amount of diagnostic output. loglevel - /// = 0 suppresses diagnostics, and increasingly-verbose messages - /// are written as loglevel increases. The messages are written to - /// a file in HTML format for viewing in a web browser. - /// @see HTML_logs - /// - /// @ingroup equil - - int equilibrate(thermo_t& s, const char* XY, int solver, - doublereal rtol, int maxsteps, int maxiter, int loglevel) { - MultiPhase* m = 0; - ChemEquil* e = 0; - bool redo = true; - int retn = -1; - int nAttempts = 0; - int retnSub = 0; - - beginLogGroup("equilibrate", loglevel); - addLogEntry("Single-phase equilibrate function"); - { - beginLogGroup("arguments"); - addLogEntry("phase",s.id()); - addLogEntry("XY",XY); - addLogEntry("solver",solver); - addLogEntry("rtol",rtol); - addLogEntry("maxsteps",maxsteps); - addLogEntry("maxiter",maxiter); - addLogEntry("loglevel",loglevel); - endLogGroup("arguments"); - } - while (redo) { - if (solver > 0) { - m = new MultiPhase; - try { - m->addPhase(&s, 1.0); - m->init(); - nAttempts++; - (void) equilibrate(*m, XY, rtol, maxsteps, maxiter, loglevel); - redo = false; - addLogEntry("MultiPhaseEquil solver succeeded."); - delete m; - retn = nAttempts; - } - catch (CanteraError err) { - addLogEntry("MultiPhaseEquil solver failed."); - delete m; - if (nAttempts < 2) { - addLogEntry("Trying single phase ChemEquil solver."); - solver = -1; - } else { - endLogGroup("equilibrate"); - throw err; - } - } - } - else { // solver <= 0 - /* - * Call the element potential solver - */ - e = new ChemEquil; - try { - e->options.maxIterations = maxsteps; - e->options.relTolerance = rtol; - nAttempts++; - retnSub = e->equilibrate(s,XY); - if (retnSub < 0) { - addLogEntry("ChemEquil solver failed."); - if (nAttempts < 2) { - addLogEntry("Trying MultiPhaseEquil solver."); - solver = 1; - } else { - throw CanteraError("equilibrate", - "Both equilibrium solvers failed"); - } - } - retn = nAttempts; - s.setElementPotentials(e->elementPotentials()); - redo = false; - delete e; - addLogEntry("ChemEquil solver succeeded."); - } - - catch (CanteraError err) { - delete e; - addLogEntry("ChemEquil solver failed."); - // If ChemEquil fails, try the MultiPhase solver - if (solver < 0) { - addLogEntry("Trying MultiPhaseEquil solver."); - solver = 1; - } - else { - redo = false; - endLogGroup("equilibrate"); - throw err; - } - } - } - } // while (redo) - /* - * We are here only for a success - */ - endLogGroup("equilibrate"); - return retn; + else { + addLogEntry("multiphase equilibrium can be done only for TP, HP, SP, or TV"); + endLogGroup("equilibrate"); + throw CanteraError("equilibrate","unsupported option"); + return -1.0; } + } + + /* + * Set a single-phase chemical solution to chemical equilibrium. + * This is a convenience function that uses one or the other of + * the two chemical equilibrium solvers. + * + * @param s The object to set to an equilibrium state + * + * @param XY An integer specifying the two properties to be held + * constant. + * + * @param solver The equilibrium solver to use. If solver = 0, + * the ChemEquil solver will be used, and if solver = 1, the + * MultiPhaseEquil solver will be used (slower than ChemEquil, + * but more stable). If solver < 0 (default, then ChemEquil will + * be tried first, and if it fails MultiPhaseEquil will be tried. + * + * @param maxsteps The maximum number of steps to take to find + * the solution. + * + * @param maxiter For the MultiPhaseEquil solver only, this is + * the maximum number of outer temperature or pressure iterations + * to take when T and/or P is not held fixed. + * + * @param loglevel Controls amount of diagnostic output. loglevel + * = 0 suppresses diagnostics, and increasingly-verbose messages + * are written as loglevel increases. The messages are written to + * a file in HTML format for viewing in a web browser. + * @see HTML_logs + * + * @ingroup equil + */ + int equilibrate(thermo_t& s, const char* XY, int solver, + doublereal rtol, int maxsteps, int maxiter, int loglevel) { + MultiPhase* m = 0; + ChemEquil* e = 0; + bool redo = true; + int retn = -1; + int nAttempts = 0; + int retnSub = 0; + + beginLogGroup("equilibrate", loglevel); + addLogEntry("Single-phase equilibrate function"); + { + beginLogGroup("arguments"); + addLogEntry("phase",s.id()); + addLogEntry("XY",XY); + addLogEntry("solver",solver); + addLogEntry("rtol",rtol); + addLogEntry("maxsteps",maxsteps); + addLogEntry("maxiter",maxiter); + addLogEntry("loglevel",loglevel); + endLogGroup("arguments"); + } + while (redo) { + if (solver > 0) { + m = new MultiPhase; + try { + m->addPhase(&s, 1.0); + m->init(); + nAttempts++; + (void) equilibrate(*m, XY, rtol, maxsteps, maxiter, loglevel); + redo = false; + addLogEntry("MultiPhaseEquil solver succeeded."); + delete m; + retn = nAttempts; + } + catch (CanteraError err) { + addLogEntry("MultiPhaseEquil solver failed."); + delete m; + if (nAttempts < 2) { + addLogEntry("Trying single phase ChemEquil solver."); + solver = -1; + } else { + endLogGroup("equilibrate"); + throw err; + } + } + } + else { // solver <= 0 + /* + * Call the element potential solver + */ + e = new ChemEquil; + try { + e->options.maxIterations = maxsteps; + e->options.relTolerance = rtol; + nAttempts++; + retnSub = e->equilibrate(s,XY); + if (retnSub < 0) { + addLogEntry("ChemEquil solver failed."); + if (nAttempts < 2) { + addLogEntry("Trying MultiPhaseEquil solver."); + solver = 1; + } else { + throw CanteraError("equilibrate", + "Both equilibrium solvers failed"); + } + } + retn = nAttempts; + s.setElementPotentials(e->elementPotentials()); + redo = false; + delete e; + addLogEntry("ChemEquil solver succeeded."); + } + + catch (CanteraError err) { + delete e; + addLogEntry("ChemEquil solver failed."); + // If ChemEquil fails, try the MultiPhase solver + if (solver < 0) { + addLogEntry("Trying MultiPhaseEquil solver."); + solver = 1; + } + else { + redo = false; + endLogGroup("equilibrate"); + throw err; + } + } + } + } // while (redo) + /* + * We are here only for a success + */ + endLogGroup("equilibrate"); + return retn; + } }