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

Cleaned up Doxygen docs for class ChemEquil

Browse source
This commit is contained in:
Ray Speth 2013-04-18 22:07:43 +00:00
parent d1e438d8dd
commit 651ce785cc
2 changed files with 106 additions and 149 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

126
include/cantera/equil/ChemEquil.h
View file

@ -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<size_t> 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<size_t> m_orderVectorElements;
std::vector<size_t> m_orderVectorSpecies;
};
extern int ChemEquil_print_lvl;

129
src/equil/ChemEquil.cpp
View file

@ -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