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

Clean up comments in other equilibrium-related classes

Browse source
This commit is contained in:
Ray Speth 2015-10-21 16:12:54 -04:00
parent 618c485c55
commit 5bf3b0fc67
6 changed files with 434 additions and 628 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -48,10 +48,9 @@ public:
int propertyPair;
/**
* Continuation flag. Set true if the calculation should be
* initialized from the last calculation. Otherwise, the
* calculation will be started from scratch and the initial
* composition and element potentials estimated.
* Continuation flag. Set true if the calculation should be initialized from
* the last calculation. Otherwise, the calculation will be started from
* scratch and the initial composition and element potentials estimated.
*/
bool contin;
};
@ -64,26 +63,22 @@ class PropertyCalculator;
*/
/**
* Class ChemEquil implements a chemical equilibrium solver for
* single-phase solutions. It is a "non-stoichiometric" solver in
* the terminology of Smith and Missen, meaning that every
* intermediate state is a valid chemical equilibrium state, but
* does not necessarily satisfy the element constraints. In
* contrast, the solver implemented in class MultiPhaseEquil uses
* a "stoichiometric" algorithm, in which each intermediate state
* satisfies the element constraints but is not a state of
* chemical equilibrium. Non-stoichiometric methods are faster
* when they converge, but stoichiometric ones tend to be more
* robust and can be used also for problems with multiple
* condensed phases. As expected, the ChemEquil solver is faster
* than MultiPhaseEquil for many single-phase equilibrium
* problems (particularly if there are only a few elements but
* very many species), but can be less stable. Problem
* situations include low temperatures where only a few species
* have non-zero mole fractions, precisely stoichiometric
* compositions (e.g. 2 H2 + O2). In general, if speed is
* important, this solver should be tried first, and if it fails
* then use MultiPhaseEquil.
* Class ChemEquil implements a chemical equilibrium solver for single-phase
* solutions. It is a "non-stoichiometric" solver in the terminology of Smith
* and Missen, meaning that every intermediate state is a valid chemical
* equilibrium state, but does not necessarily satisfy the element constraints.
* In contrast, the solver implemented in class MultiPhaseEquil uses a
* "stoichiometric" algorithm, in which each intermediate state satisfies the
* element constraints but is not a state of chemical equilibrium. Non-
* stoichiometric methods are faster when they converge, but stoichiometric ones
* tend to be more robust and can be used also for problems with multiple
* condensed phases. As expected, the ChemEquil solver is faster than
* MultiPhaseEquil for many single-phase equilibrium problems (particularly if
* there are only a few elements but very many species), but can be less stable.
* Problem situations include low temperatures where only a few species have
* non-zero mole fractions, precisely stoichiometric compositions (e.g. 2 H2 +
* O2). In general, if speed is important, this solver should be tried first,
* and if it fails then use MultiPhaseEquil.
* @ingroup equil
*/
class ChemEquil
@ -93,32 +88,33 @@ public:
//! 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.
* 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(thermo_t& s);
virtual ~ChemEquil();
/*!
* Equilibrate a phase, holding the elemental composition fixed
* at the initial value found within the ThermoPhase object *s*.
* 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.
* The value of two 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.
* Compute the equilibrium composition for two 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.
* The two specified properties are obtained by querying the ThermoPhase
* object. The properties must be already contained within the current
* thermodynamic state of the system.
*
* @param s phase object to be equilibrated
* @param XY property pair to hold constant
@ -146,10 +142,9 @@ public:
protected:
//! Pointer to the ThermoPhase object used to initialize this object.
/*!
* This ThermoPhase object must be compatible with the ThermoPhase
* objects input from the equilibrate function. Currently, this
* means that the 2 ThermoPhases have to have consist of the same
* species and elements.
* This ThermoPhase object must be compatible with the ThermoPhase objects
* input from the equilibrate function. Currently, this means that the 2
* ThermoPhases have to have consist of the same species and elements.
*/
thermo_t* m_phase;
@ -159,14 +154,14 @@ protected:
}
/*!
* Prepare for equilibrium calculations.
* @param s object representing the solution phase.
* 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.
* Set mixture to an equilibrium state consistent with specified element
* potentials and temperature.
*
* @param s mixture to be updated
* @param x vector of non-dimensional element potentials
@ -267,51 +262,36 @@ protected:
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.
*/
//! This is equal to the rank of the stoichiometric coefficient matrix when
//! it is computed. It's initialized to #m_mm.
size_t m_nComponents;
std::unique_ptr<PropertyCalculator<thermo_t> > m_p1, m_p2;
/**
* Current value of the mole fractions in the single phase.
* -> length = #m_kk.
*/
//! Current value of the mole fractions in the single phase. length = #m_kk.
vector_fp m_molefractions;
/**
* Current value of the dimensional element potentials
* -> length = #m_mm
*/
//! Current value of the dimensional element potentials. length = #m_mm
vector_fp m_lambda;
/*
* Current value of the sum of the element abundances given the
* current element potentials.
*/
//! Current value of the sum of the element abundances given the current
//! element potentials.
doublereal m_elementTotalSum;
/*
* Current value of the element mole fractions. Note these aren't
* the goal element mole fractions.
*/
//! Current value of the element mole fractions. Note these aren't the goal
//! element mole fractions.
vector_fp m_elementmolefracs;
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];
*/
//! Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];
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.
*/
//! Index of the element id corresponding to the electric charge of each
//! species. Equal to -1 if there is no such element id.
size_t m_eloc;
vector_fp m_startSoln;
@ -319,11 +299,9 @@ 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
* pressure of the solution (the star standard state).
*/
//! Dimensionless values of the Gibbs free energy for the standard state of
//! each species, at the temperature and pressure of the solution (the star
//! standard state).
vector_fp m_muSS_RT;
std::vector<size_t> m_component;

153
include/cantera/equil/MultiPhase.h
View file

@ -145,7 +145,7 @@ public:
/*!
* @param kGlob global species index
* @param mGlob global element index
* @return returns the number of atoms.
* @returns the number of atoms.
*/
doublereal nAtoms(const size_t kGlob, const size_t mGlob) const;
@ -176,9 +176,7 @@ public:
//! Returns the index, given the phase name
/*!
* @param pName Name of the phase
*
* @return returns the index. A value of -1 means
* the phase isn't in the object.
* @returns the index. A value of -1 means the phase isn't in the object.
*/
int phaseIndex(const std::string& pName) const;
@ -236,9 +234,10 @@ public:
* @param speciesName Species Name
* @param phaseName Phase Name
*
* @return returns the global index
* If the species or phase name is not recognized, this routine throws
* a CanteraError.
* @returns the global index
*
* If the species or phase name is not recognized, this routine throws a
* CanteraError.
*/
size_t speciesIndex(const std::string& speciesName, const std::string& phaseName);
@ -273,14 +272,14 @@ public:
*/
doublereal elementMoles(size_t m) const;
//! Returns a vector of Chemical potentials.
//! Returns a vector of Chemical potentials.
/*!
* Write into array \a mu the chemical potentials of all species
* [J/kmol]. The chemical potentials are related to the activities by
* Write into array \a mu the chemical potentials of all species
* [J/kmol]. The chemical potentials are related to the activities by
*
* \f$
* \f$
* \mu_k = \mu_k^0(T, P) + RT \ln a_k.
* \f$.
* \f$.
*
* @param mu Chemical potential vector. Length = num global species. Units
* = J/kmol.
@ -372,17 +371,17 @@ public:
//! Set the state of the underlying ThermoPhase objects in one call
/*!
* @param T Temperature of the system (kelvin)
* @param Pres pressure of the system (pascal)
* @param T Temperature of the system (kelvin)
* @param Pres pressure of the system (pascal)
*/
void setState_TP(const doublereal T, const doublereal Pres);
//! Set the state of the underlying ThermoPhase objects in one call
/*!
* @param T Temperature of the system (kelvin)
* @param Pres pressure of the system (pascal)
* @param Moles Vector of mole numbers of all the species in all the phases
* (kmol)
* @param T Temperature of the system (kelvin)
* @param Pres pressure of the system (pascal)
* @param Moles Vector of mole numbers of all the species in all the phases
* (kmol)
*/
void setState_TPMoles(const doublereal T, const doublereal Pres, const doublereal* Moles);
@ -393,8 +392,8 @@ public:
/// The total mixture volume [m^3].
/*!
* Returns the cumulative sum of the volumes of all the
* phases in the mixture.
* Returns the cumulative sum of the volumes of all the phases in the
* mixture.
*/
doublereal volume() const;
@ -428,8 +427,8 @@ public:
return m_np;
}
//! Return true is species \a kGlob is a species in a
//! multicomponent solution phase.
//! Return true is species \a kGlob is a species in a multicomponent
//! solution phase.
/*!
* @param kGlob index of the global species
*/
@ -438,7 +437,7 @@ public:
//! Returns the phase index of the Kth "global" species
/*!
* @param kGlob Global species index.
* @return Returns the index of the owning phase.
* @returns the index of the owning phase.
*/
size_t speciesPhaseIndex(const size_t kGlob) const;
@ -495,8 +494,8 @@ public:
//! Adds moles of a certain species to the mixture
/*!
* @param indexS Index of the species in the MultiPhase object
* @param addedMoles Value of the moles that are added to the species.
* @param indexS Index of the species in the MultiPhase object
* @param addedMoles Value of the moles that are added to the species.
*/
void addSpeciesMoles(const int indexS, const doublereal addedMoles);
@ -508,8 +507,8 @@ public:
*/
void getElemAbundances(doublereal* elemAbundances) const;
//! Return true if the phase \a p has valid thermo data for
//! the current temperature.
//! Return true if the phase \a p has valid thermo data for the current
//! temperature.
/*!
* @param p Index of the phase.
*/
@ -517,31 +516,30 @@ public:
// These methods are meant for internal use.
//! Update the locally-stored composition within this object
//! to match the current compositions of the phase objects.
//! Update the locally-stored composition within this object to match the
//! current compositions of the phase objects.
/*!
* Query the underlying ThermoPhase objects for their mole fractions and
* fill in the mole fraction vector of this current object. Adjust
* element compositions within this object to match.
* Query the underlying ThermoPhase objects for their mole fractions and
* fill in the mole fraction vector of this current object. Adjust element
* compositions within this object to match.
*
* This is an upload operation in the sense that we are taking downstream
* information (ThermoPhase object info) and applying it to an upstream
* object (MultiPhase object).
* This is an upload operation in the sense that we are taking downstream
* information (ThermoPhase object info) and applying it to an upstream
* object (MultiPhase object).
*/
void uploadMoleFractionsFromPhases();
//! Set the states of the phase objects to the locally-stored
//! state within this MultiPhase object.
/*!
* This method sets each phase to the mixture temperature and pressure,
* and sets the phase mole fractions based on the mixture mole numbers.
* This method sets each phase to the mixture temperature and pressure,
* and sets the phase mole fractions based on the mixture mole numbers.
*
* This is an download operation in the sense that we are taking
* upstream object information (MultiPhase object) and
* applying it to downstream objects (ThermoPhase object information)
* This is an download operation in the sense that we are taking upstream
* object information (MultiPhase object) and applying it to downstream
* objects (ThermoPhase object information)
*
* Therefore, the term, "update", is appropriate for a downstream
* operation.
* Therefore, the term, "update", is appropriate for a downstream operation.
*/
void updatePhases() const;
@ -552,15 +550,14 @@ private:
//! Set the mixture to a state of chemical equilibrium using the
//! MultiPhaseEquil solver.
/*!
* @param XY Integer flag specifying properties to hold fixed.
* @param err Error tolerance for \f$\Delta \mu/RT \f$ for all reactions.
* Also used as the relative error tolerance for the outer
* loop.
* @param maxsteps Maximum number of steps to take in solving the fixed
* TP problem.
* @param maxiter Maximum number of "outer" iterations for problems holding
* fixed something other than (T,P).
* @param loglevel Level of diagnostic output
* @param XY Integer flag specifying properties to hold fixed.
* @param err Error tolerance for \f$\Delta \mu/RT \f$ for all reactions.
* Also used as the relative error tolerance for the outer loop.
* @param maxsteps Maximum number of steps to take in solving the fixed TP
* problem.
* @param maxiter Maximum number of "outer" iterations for problems holding
* fixed something other than (T,P).
* @param loglevel Level of diagnostic output
*/
double equilibrate_MultiPhaseEquil(int XY, doublereal err, int maxsteps,
int maxiter, int loglevel);
@ -576,16 +573,14 @@ private:
//! Global Stoichiometric Coefficient array
/*!
* This is a two dimensional array m_atoms(m, k). The first index is the
* global element index. The second index, k, is the global species
* index. The value is the number of atoms of type m in species k.
* This is a two dimensional array m_atoms(m, k). The first index is the
* global element index. The second index, k, is the global species index.
* The value is the number of atoms of type m in species k.
*/
DenseMatrix m_atoms;
/**
* Locally stored vector of mole fractions of all species
* comprising the MultiPhase object.
*/
//! Locally stored vector of mole fractions of all species comprising the
//! MultiPhase object.
vector_fp m_moleFractions;
//! Mapping between the global species number and the phase ID
@ -603,23 +598,15 @@ private:
*/
std::vector<size_t> m_spstart;
//! String names of the global elements
/*!
* This has a length equal to the number of global elements.
*/
//! String names of the global elements. This has a length equal to the
//! number of global elements.
std::vector<std::string> m_enames;
//! Atomic number of each element
/*!
* This is the atomic number of each global element.
*/
//! Atomic number of each global element.
vector_int m_atomicNumber;
//! Vector of species names in the problem
/*!
* Vector is over all species defined in the object,
* the global species index.
*/
//! Vector of species names in the problem. Vector is over all species
//! defined in the object, the global species index.
std::vector<std::string> m_snames;
//! Returns the global element index, given the element string name
@ -646,10 +633,8 @@ private:
//! True if the init() routine has been called, and the MultiPhase frozen
bool m_init;
//! Global ID of the element corresponding to the electronic charge.
/*!
* If there is none, then this is equal to -1
*/
//! Global ID of the element corresponding to the electronic charge. If
//! there is none, then this is equal to -1
size_t m_eloc;
//! Vector of bools indicating whether temperatures are ok for phases.
@ -659,17 +644,12 @@ private:
*/
mutable std::vector<bool> m_temp_OK;
//! Minimum temperature for which thermo parameterizations are valid
/*!
* Stoichiometric phases are ignored in this determination.
* units Kelvin
*/
//! Minimum temperature for which thermo parameterizations are valid.
//! Stoichiometric phases are ignored in this determination. units Kelvin
doublereal m_Tmin;
//! Minimum temperature for which thermo parameterizations are valid
/*!
* Stoichiometric phases are ignored in this determination. units Kelvin
*/
//! Minimum temperature for which thermo parameterizations are valid.
//! Stoichiometric phases are ignored in this determination. units Kelvin
doublereal m_Tmax;
//! Vector of element abundances
@ -711,9 +691,8 @@ inline std::ostream& operator<<(std::ostream& s, MultiPhase& x)
* currently a linear combination of the previous components. Then, calculate
* the stoichiometric coefficient matrix for that basis.
*
* Calculates the identity of the component species in the mechanism.
* Rearranges the solution data to put the component data at the
* front of the species list.
* Calculates the identity of the component species in the mechanism. Rearranges
* the solution data to put the component data at the front of the species list.
*
* Then, calculates SC(J,I) the formation reactions for all noncomponent
* species in the mechanism.

296
src/equil/BasisOptimize.cpp
View file

@ -31,20 +31,16 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
size_t j, jj, k=0, kk, l, i, jl, ml;
std::string ename;
std::string sname;
/*
* Get the total number of elements defined in the multiphase object
*/
// Get the total number of elements defined in the multiphase object
size_t ne = mphase->nElements();
/*
* Get the total number of species in the multiphase object
*/
// Get the total number of species in the multiphase object
size_t nspecies = mphase->nSpecies();
doublereal tmp;
doublereal const USEDBEFORE = -1;
/*
* Perhaps, initialize the element ordering
*/
// Perhaps, initialize the element ordering
if (orderVectorElements.size() < ne) {
orderVectorElements.resize(ne);
for (j = 0; j < ne; j++) {
@ -52,9 +48,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
}
}
/*
* Perhaps, initialize the species ordering
*/
// Perhaps, initialize the species ordering
if (orderVectorSpecies.size() != nspecies) {
orderVectorSpecies.resize(nspecies);
for (k = 0; k < nspecies; k++) {
@ -101,27 +95,20 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
}
}
/*
* Calculate the maximum value of the number of components possible
* It's equal to the minimum of the number of elements and the
* number of total species.
*/
// Calculate the maximum value of the number of components possible. It's
// equal to the minimum of the number of elements and the number of total
// species.
size_t nComponents = std::min(ne, nspecies);
size_t nNonComponents = nspecies - nComponents;
/*
* Set this return variable to false
*/
// Set this return variable to false
*usedZeroedSpecies = false;
/*
* Create an array of mole numbers
*/
// Create an array of mole numbers
vector_fp molNum(nspecies,0.0);
mphase->getMoles(molNum.data());
/*
* Other workspace
*/
// Other workspace
vector_fp sm(ne*ne, 0.0);
vector_fp ss(ne, 0.0);
vector_fp sa(ne, 0.0);
@ -129,9 +116,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
formRxnMatrix.resize(nspecies*ne, 0.0);
}
/*
* For debugging purposes keep an unmodified copy of the array.
*/
// For debugging purposes keep an unmodified copy of the array.
vector_fp molNumBase;
if (DEBUG_MODE_ENABLED) {
molNumBase = molNum;
@ -139,19 +124,16 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
double molSave = 0.0;
size_t jr = 0;
/*
* Top of a loop of some sort based on the index JR. JR is the
* current number of component species found.
*/
// Top of a loop of some sort based on the index JR. JR is the current
// number of component species found.
while (jr < nComponents) {
// Top of another loop point based on finding a linearly independent
// species
while (true) {
/*
* Search the remaining part of the mole number vector, molNum
* for the largest remaining species. Return its identity.
* kk is the raw number. k is the orderVectorSpecies index.
*/
// Search the remaining part of the mole number vector, molNum for
// the largest remaining species. Return its identity. kk is the raw
// number. k is the orderVectorSpecies index.
kk = max_element(molNum.begin(), molNum.end()) - molNum.begin();
for (j = 0; j < nspecies; j++) {
if (orderVectorSpecies[j] == kk) {
@ -166,40 +148,33 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
if (molNum[kk] == 0.0) {
*usedZeroedSpecies = true;
}
/*
* If the largest molNum is negative, then we are done.
*/
// If the largest molNum is negative, then we are done.
if (molNum[kk] == USEDBEFORE) {
nComponents = jr;
nNonComponents = nspecies - nComponents;
break;
}
/*
* Assign a small negative number to the component that we have
* just found, in order to take it out of further consideration.
*/
// Assign a small negative number to the component that we have
// just found, in order to take it out of further consideration.
#ifdef DEBUG_MODE
molSave = molNum[kk];
#endif
molNum[kk] = USEDBEFORE;
// **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ******
/*
* Modified Gram-Schmidt Method, p. 202 Dalquist
* QR factorization of a matrix without row pivoting.
*/
// CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES
// Modified Gram-Schmidt Method, p. 202 Dalquist
// QR factorization of a matrix without row pivoting.
jl = jr;
for (j = 0; j < ne; ++j) {
jj = orderVectorElements[j];
sm[j + jr*ne] = mphase->nAtoms(kk,jj);
}
if (jl > 0) {
/*
* Compute the coefficients of JA column of the
* the upper triangular R matrix, SS(J) = R_J_JR
* (this is slightly different than Dalquist)
* R_JA_JA = 1
*/
// Compute the coefficients of JA column of the the upper
// triangular R matrix, SS(J) = R_J_JR (this is slightly
// different than Dalquist) R_JA_JA = 1
for (j = 0; j < jl; ++j) {
ss[j] = 0.0;
for (i = 0; i < ne; ++i) {
@ -207,31 +182,31 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
}
ss[j] /= sa[j];
}
/*
* Now make the new column, (*,JR), orthogonal to the
* previous columns
*/
// Now make the new column, (*,JR), orthogonal to the previous
// columns
for (j = 0; j < jl; ++j) {
for (l = 0; l < ne; ++l) {
sm[l + jr*ne] -= ss[j] * sm[l + j*ne];
}
}
}
/*
* Find the new length of the new column in Q.
* It will be used in the denominator in future row calcs.
*/
// Find the new length of the new column in Q.
// It will be used in the denominator in future row calcs.
sa[jr] = 0.0;
for (ml = 0; ml < ne; ++ml) {
tmp = sm[ml + jr*ne];
sa[jr] += tmp * tmp;
}
// **** IF NORM OF NEW ROW .LT. 1E-3 REJECT **********
// IF NORM OF NEW ROW .LT. 1E-3 REJECT
if (sa[jr] > 1.0e-6) {
break;
}
}
// **** REARRANGE THE DATA ******************
// REARRANGE THE DATA
if (jr != k) {
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) {
kk = orderVectorSpecies[k];
@ -252,38 +227,36 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
return nComponents;
}
// **** EVALUATE THE STOICHIOMETRY **********************
/*
* Formulate the matrix problem for the stoichiometric
* coefficients. CX + B = 0
* C will be an nc x nc matrix made up of the formula
* vectors for the components. Each component's formula
* vector is a column. The rows are the elements.
* n RHS's will be solved for. Thus, B is an nc x n
* matrix.
*
* BIG PROBLEM 1/21/99:
*
* This algorithm makes the assumption that the
* first nc rows of the formula matrix aren't rank deficient.
* However, this might not be the case. For example, assume
* that the first element in FormulaMatrix[] is argon. Assume that
* no species in the matrix problem actually includes argon.
* Then, the first row in sm[], below will be identically
* zero. bleh.
* What needs to be done is to perform a rearrangement
* of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and gai, such
* that the first nc elements form in combination with the
* nc components create an invertible sm[]. not a small
* project, but very doable.
* An alternative would be to turn the matrix problem
* below into an ne x nc problem, and do QR elimination instead
* of Gauss-Jordan elimination.
* Note the rearrangement of elements need only be done once
* in the problem. It's actually very similar to the top of
* this program with ne being the species and nc being the
* elements!!
*/
// EVALUATE THE STOICHIOMETRY
//
// Formulate the matrix problem for the stoichiometric
// coefficients. CX + B = 0
//
// C will be an nc x nc matrix made up of the formula vectors for the
// components. Each component's formula vector is a column. The rows are the
// elements.
//
// n RHS's will be solved for. Thus, B is an nc x n matrix.
//
// BIG PROBLEM 1/21/99:
//
// This algorithm makes the assumption that the first nc rows of the formula
// matrix aren't rank deficient. However, this might not be the case. For
// example, assume that the first element in FormulaMatrix[] is argon.
// Assume that no species in the matrix problem actually includes argon.
// Then, the first row in sm[], below will be identically zero. bleh.
//
// What needs to be done is to perform a rearrangement of the ELEMENTS ->
// i.e. rearrange, FormulaMatrix, sp, and gai, such that the first nc
// elements form in combination with the nc components create an invertible
// sm[]. not a small project, but very doable.
//
// An alternative would be to turn the matrix problem below into an ne x nc
// problem, and do QR elimination instead of Gauss-Jordan elimination.
//
// Note the rearrangement of elements need only be done once in the problem.
// It's actually very similar to the top of this program with ne being the
// species and nc being the elements!!
for (k = 0; k < nComponents; ++k) {
kk = orderVectorSpecies[k];
for (j = 0; j < nComponents; ++j) {
@ -337,9 +310,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
writelogf(" --- %3d (%3d) ", k, kk);
writelogf("%-10.10s", mphase->speciesName(kk));
writelogf("|%10.3g|", molNumBase[kk]);
/*
* Print the negative of formRxnMatrix[]; it's easier to interpret.
*/
// Print the negative of formRxnMatrix[]; it's easier to interpret.
for (j = 0; j < nComponents; j++) {
writelogf(" %6.2f", - formRxnMatrix[j + i * ne]);
}
@ -355,21 +326,17 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
return nComponents;
} // basopt()
/**
* Print a string within a given space limit. This routine limits the amount of
* the string that will be printed to a maximum of "space" characters.
*
* str = String -> must be null terminated.
* space = space limit for the printing.
* alignment = 0 centered
* 1 right aligned
* 2 left aligned
*/
static void print_stringTrunc(const char* str, int space, int alignment)
/***********************************************************************
* vcs_print_stringTrunc():
*
* Print a string within a given space limit. This routine
* limits the amount of the string that will be printed to a
* maximum of "space" characters.
*
* str = String -> must be null terminated.
* space = space limit for the printing.
* alignment = 0 centered
* 1 right aligned
* 2 left aligned
***********************************************************************/
{
int i, ls=0, rs=0;
int len = static_cast<int>(strlen(str));
@ -408,9 +375,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
size_t j, k, l, i, jl, ml, jr, ielem, jj, kk=0;
size_t nelements = mphase->nElements();
std::string ename;
/*
* Get the total number of species in the multiphase object
*/
// Get the total number of species in the multiphase object
size_t nspecies = mphase->nSpecies();
double test = -1.0E10;
@ -425,9 +390,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
writelog(" --- and to rearrange the element ordering once\n");
}
/*
* Perhaps, initialize the element ordering
*/
// Perhaps, initialize the element ordering
if (orderVectorElements.size() < nelements) {
orderVectorElements.resize(nelements);
for (j = 0; j < nelements; j++) {
@ -435,11 +398,8 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
}
}
/*
* Perhaps, initialize the species ordering. However, this is
* dangerous, as this ordering is assumed to yield the
* component species for the problem
*/
// Perhaps, initialize the species ordering. However, this is dangerous, as
// this ordering is assumed to yield the component species for the problem
if (orderVectorSpecies.size() != nspecies) {
orderVectorSpecies.resize(nspecies);
for (k = 0; k < nspecies; k++) {
@ -447,12 +407,9 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
}
}
/*
* If the elementAbundances aren't input, just create a fake one
* based on summing the column of the stoich matrix.
* This will force elements with zero species to the
* end of the element ordering.
*/
// If the elementAbundances aren't input, just create a fake one based on
// summing the column of the stoich matrix. This will force elements with
// zero species to the end of the element ordering.
vector_fp eAbund(nelements,0.0);
if (elementAbundances.size() != nelements) {
for (j = 0; j < nelements; j++) {
@ -470,23 +427,16 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
vector_fp ss(nelements,0.0);
vector_fp sm(nelements*nelements,0.0);
/*
* Top of a loop of some sort based on the index JR. JR is the
* current number independent elements found.
*/
// Top of a loop of some sort based on the index JR. JR is the current
// number independent elements found.
jr = 0;
while (jr < nComponents) {
/*
* Top of another loop point based on finding a linearly
* independent element
*/
// Top of another loop point based on finding a linearly independent
// element
while (true) {
/*
* Search the element vector. We first locate elements that
* are present in any amount. Then, we locate elements that
* are not present in any amount.
* Return its identity in K.
*/
// Search the element vector. We first locate elements that are
// present in any amount. Then, we locate elements that are not
// present in any amount. Return its identity in K.
k = nelements;
for (ielem = jr; ielem < nelements; ielem++) {
kk = orderVectorElements[ielem];
@ -512,37 +462,30 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
throw CanteraError("ElemRearrange", "Required number of elements not found.");
}
/*
* Assign a large negative number to the element that we have
* just found, in order to take it out of further consideration.
*/
// Assign a large negative number to the element that we have
// just found, in order to take it out of further consideration.
eAbund[kk] = test;
// **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX
// **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ******
/*
* Modified Gram-Schmidt Method, p. 202 Dalquist
* QR factorization of a matrix without row pivoting.
*/
// CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX
// LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX
// Modified Gram-Schmidt Method, p. 202 Dalquist
// QR factorization of a matrix without row pivoting.
jl = jr;
/*
* Fill in the row for the current element, k, under consideration
* The row will contain the Formula matrix value for that element
* with respect to the vector of component species.
* (note j and k indices are flipped compared to the previous routine)
*/
// Fill in the row for the current element, k, under consideration
// The row will contain the Formula matrix value for that element
// with respect to the vector of component species. (note j and k
// indices are flipped compared to the previous routine)
for (j = 0; j < nComponents; ++j) {
jj = orderVectorSpecies[j];
kk = orderVectorElements[k];
sm[j + jr*nComponents] = mphase->nAtoms(jj,kk);
}
if (jl > 0) {
/*
* Compute the coefficients of JA column of the
* the upper triangular R matrix, SS(J) = R_J_JR
* (this is slightly different than Dalquist)
* R_JA_JA = 1
*/
// Compute the coefficients of JA column of the the upper
// triangular R matrix, SS(J) = R_J_JR (this is slightly
// different than Dalquist) R_JA_JA = 1
for (j = 0; j < jl; ++j) {
ss[j] = 0.0;
for (i = 0; i < nComponents; ++i) {
@ -550,10 +493,9 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
}
ss[j] /= sa[j];
}
/*
* Now make the new column, (*,JR), orthogonal to the
* previous columns
*/
// Now make the new column, (*,JR), orthogonal to the
// previous columns
for (j = 0; j < jl; ++j) {
for (l = 0; l < nComponents; ++l) {
sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents];
@ -561,21 +503,19 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
}
}
/*
* Find the new length of the new column in Q.
* It will be used in the denominator in future row calcs.
*/
// Find the new length of the new column in Q.
// It will be used in the denominator in future row calcs.
sa[jr] = 0.0;
for (ml = 0; ml < nComponents; ++ml) {
double tmp = sm[ml + jr*nComponents];
sa[jr] += tmp * tmp;
}
// **** IF NORM OF NEW ROW .LT. 1E-6 REJECT **********
// IF NORM OF NEW ROW .LT. 1E-6 REJECT
if (sa[jr] > 1.0e-6) {
break;
}
}
// **** REARRANGE THE DATA ******************
// REARRANGE THE DATA
if (jr != k) {
if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
kk = orderVectorElements[k];

344
src/equil/ChemEquil.cpp
View file

@ -201,11 +201,9 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal,
for (size_t m = 0; m < m_nComponents; m++) {
m_component[m] = e.componentIndex(m);
}
/*
* Update the current values of the temp, density, and
* mole fraction, and element abundance vectors kept
* within the ChemEquil object.
*/
// Update the current values of the temp, density, and mole fraction,
// and element abundance vectors kept within the ChemEquil object.
update(s);
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
@ -347,9 +345,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
vector_fp state;
s.saveState(state);
/*
* Check Compatibility
*/
// Check Compatibility
if (m_mm != s.nElements() || m_kk != s.nSpecies()) {
throw CanteraError("ChemEquil::equilibrate ERROR",
"Input ThermoPhase is incompatible with initialization");
@ -406,11 +402,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
}
}
/*
* Before we do anything to change the ThermoPhase object,
* we calculate and store the two specified thermodynamic
* properties that we are after.
*/
// Before we do anything to change the ThermoPhase object, we calculate and
// store the two specified thermodynamic properties that we are after.
xval = m_p1->value(s);
yval = m_p2->value(s);
@ -420,13 +413,10 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
vector_fp x(nvar, -102.0); // solution vector
vector_fp res_trial(nvar, 0.0); // residual
/*
* Replace one of the element abundance fraction equations
* with the specified property calculation.
*
* We choose the equation of the element with the highest element
* abundance.
*/
// Replace one of the element abundance fraction equations with the
// specified property calculation.
//
// We choose the equation of the element with the highest element abundance.
size_t m;
tmp = -1.0;
for (size_t im = 0; im < m_nComponents; im++) {
@ -441,20 +431,17 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
"Element Abundance Vector is zeroed");
}
// start with a composition with everything non-zero. Note
// that since we have already save the target element moles,
// changing the composition at this point only affects the
// starting point, not the final solution.
// start with a composition with everything non-zero. Note that since we
// have already save the target element moles, changing the composition at
// this point only affects the starting point, not the final solution.
vector_fp xmm(m_kk, 0.0);
for (size_t k = 0; k < m_kk; k++) {
xmm[k] = s.moleFraction(k) + 1.0E-32;
}
s.setMoleFractions(xmm.data());
/*
* Update the internally stored values of m_temp,
* m_dens, and the element mole fractions.
*/
// Update the internally stored values of m_temp, m_dens, and the element
// mole fractions.
update(s);
doublereal tmaxPhase = s.maxTemp();
@ -472,10 +459,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
doublereal slope, phigh, plow, pval, dt;
// first get the property values at the upper and lower
// temperature limits. Since p1 (h, s, or u) is monotonic
// in T, these values determine the upper and lower
// bounnds (phigh, plow) for p1.
// first get the property values at the upper and lower temperature
// limits. Since p1 (h, s, or u) is monotonic in T, these values
// determine the upper and lower bounnds (phigh, plow) for p1.
s.setTemperature(tmax);
setInitialMoles(s, elMolesGoal, loglevel - 1);
@ -495,10 +481,10 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
setInitialMoles(s, elMolesGoal, loglevel - 1);
pval = m_p1->value(s);
// If this value of p1 is greater than the specified
// property value, then the current temperature is too
// high. Use it as the new upper bound. Otherwise, it
// is too low, so use it as the new lower bound.
// If this value of p1 is greater than the specified property value,
// then the current temperature is too high. Use it as the new upper
// bound. Otherwise, it is too low, so use it as the new lower
// bound.
if (pval > xval) {
tmax = t0;
phigh = pval;
@ -534,11 +520,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
setInitialMoles(s, elMolesGoal,loglevel);
/*
* If requested, get the initial estimate for the
* chemical potentials from the ThermoPhase object
* itself. Or else, create our own estimate.
*/
// If requested, get the initial estimate for the chemical potentials from
// the ThermoPhase object itself. Or else, create our own estimate.
if (useThermoPhaseElementPotentials) {
bool haveEm = s.getElementPotentials(x.data());
if (haveEm) {
@ -552,44 +535,32 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
estimateElementPotentials(s, x, elMolesGoal);
}
} else {
/*
* Calculate initial estimates of the element potentials.
* This algorithm uese the MultiPhaseEquil object's
* initialization capabilities to calculate an initial
* estimate of the mole fractions for a set of linearly
* independent component species. Then, the element
* potentials are solved for based on the chemical
* potentials of the component species.
*/
// Calculate initial estimates of the element potentials. This algorithm
// uese the MultiPhaseEquil object's initialization capabilities to
// calculate an initial estimate of the mole fractions for a set of
// linearly independent component species. Then, the element potentials
// are solved for based on the chemical potentials of the component
// species.
estimateElementPotentials(s, x, elMolesGoal);
}
/*
* Do a better estimate of the element potentials.
* We have found that the current estimate may not be good
* enough to avoid drastic numerical issues associated with
* the use of a numerically generated Jacobian.
*
* The Brinkley algorithm assumes a constant T, P system
* and uses a linearized analytical Jacobian that turns out
* to be very stable.
*/
// Do a better estimate of the element potentials. We have found that the
// current estimate may not be good enough to avoid drastic numerical issues
// associated with the use of a numerically generated Jacobian.
//
// The Brinkley algorithm assumes a constant T, P system and uses a
// linearized analytical Jacobian that turns out to be very stable.
int info = estimateEP_Brinkley(s, x, elMolesGoal);
if (info == 0) {
setToEquilState(s, x, s.temperature());
}
/*
* Install the log(temp) into the last solution unknown
* slot.
*/
// Install the log(temp) into the last solution unknown slot.
x[m_mm] = log(s.temperature());
/*
* Setting the max and min values for x[]. Also, if element
* abundance vector is zero, setting x[] to -1000. This
* effectively zeroes out all species containing that element.
*/
// Setting the max and min values for x[]. Also, if element abundance vector
// is zero, setting x[] to -1000. This effectively zeroes out all species
// containing that element.
vector_fp above(nvar);
vector_fp below(nvar);
for (m = 0; m < mm; m++) {
@ -599,10 +570,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
x[m] = -1000.0;
}
}
/*
* Set the temperature bounds to be 25 degrees different than the max and min
* temperatures.
*/
// Set the temperature bounds to be 25 degrees different than the max and
// min temperatures.
above[mm] = log(s.maxTemp() + 25.0);
below[mm] = log(s.minTemp() - 25.0);
@ -625,19 +595,17 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
for (m = 0; m < nvar; m++) {
double tval = options.relTolerance;
if (m < mm) {
/*
* Special case convergence requirements for electron element.
* This is a special case because the element coefficients may
* be both positive and negative. And, typically they sum to 0.0.
* Therefore, there is no natural absolute value for this quantity.
* We supply the absolute value tolerance here. Note, this is
* made easier since the element abundances are normalized to one
* within this routine.
*
* Note, the 1.0E-13 value was recently relaxed from 1.0E-15, because
* convergence failures were found to occur for the lower value
* at small pressure (0.01 pascal).
*/
// Special case convergence requirements for electron element.
// This is a special case because the element coefficients may
// be both positive and negative. And, typically they sum to
// 0.0. Therefore, there is no natural absolute value for this
// quantity. We supply the absolute value tolerance here. Note,
// this is made easier since the element abundances are
// normalized to one within this routine.
//
// Note, the 1.0E-13 value was recently relaxed from 1.0E-15,
// because convergence failures were found to occur for the
// lower value at small pressure (0.01 pascal).
if (m == m_eloc) {
tval = elMolesGoal[m] * options.relTolerance + options.absElemTol
+ 1.0E-13;
@ -659,10 +627,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
if (m_eloc != npos) {
adjustEloc(s, elMolesGoal);
}
/*
* Save the calculated and converged element potentials
* to the original ThermoPhase object.
*/
// Save the calculated and converged element potentials to the
// original ThermoPhase object.
s.setElementPotentials(m_lambda);
if (s.temperature() > s.maxTemp() + 1.0 ||
s.temperature() < s.minTemp() - 1.0) {
@ -705,9 +672,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
oldf = f;
scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0);
/*
* Solve the system
*/
// Solve the system
try {
info = solve(jac, res_trial.data());
} catch (CanteraError& err) {
@ -776,9 +741,8 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx,
{
double damp;
/*
* Carry out a delta damping approach on the dimensionless element potentials.
*/
// Carry out a delta damping approach on the dimensionless element
// potentials.
damp = 1.0;
for (size_t m = 0; m < m_mm; m++) {
if (m == m_eloc) {
@ -798,9 +762,7 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx,
}
}
/*
* Update the solution unknown
*/
// Update the solution unknown
for (size_t m = 0; m < x.size(); m++) {
x[m] = oldx[m] + damp * step[m];
}
@ -832,11 +794,8 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x,
} else if (n >= m_nComponents) {
resid[m] = x[m];
} else {
/*
* Change the calculation for small element number, using
* L'Hopital's rule.
* The log formulation is unstable.
*/
// Change the calculation for small element number, using
// L'Hopital's rule. The log formulation is unstable.
if (elmFracGoal[m] < 1.0E-10 || elmFrac[m] < 1.0E-10 || m == m_eloc) {
resid[m] = elmFracGoal[m] - elmFrac[m];
} else {
@ -907,10 +866,9 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t,
{
double n_t_calc = 0.0;
double tmp;
/*
* Calculate the activity coefficients of the solution, at the
* previous solution state.
*/
// Calculate the activity coefficients of the solution, at the previous
// solution state.
vector_fp actCoeff(m_kk, 1.0);
s.setMoleFractions(Xmol_i_calc.data());
s.setPressure(pressureConst);
@ -941,11 +899,8 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t,
int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
vector_fp& elMoles)
{
/*
* Before we do anything, we will save the state of the solution.
* Then, if things go drastically wrong, we will restore the
* saved state.
*/
// Before we do anything, we will save the state of the solution. Then, if
// things go drastically wrong, we will restore the saved state.
vector_fp state;
s.saveState(state);
double tmp, sum;
@ -969,11 +924,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
vector_fp resid(m_mm+1, 0.0);
vector_int lumpSum(m_mm+1, 0);
/*
* Get the nondimensional Gibbs functions for the species
* at their standard states of solution at the current T and P
* of the solution.
*/
// Get the nondimensional Gibbs functions for the species at their standard
// states of solution at the current T and P of the solution.
s.getGibbs_RT(m_muSS_RT.data());
vector_fp eMolesCalc(m_mm, 0.0);
@ -1042,27 +994,18 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
/*
* -------------------------------------------------------------------
* Main Loop.
*/
// Main Loop.
for (int iter = 0; iter < 20* options.maxIterations; iter++) {
/*
* Save the old solution
*/
// Save the old solution
for (m = 0; m < m_mm; m++) {
x_old[m] = x[m];
}
x_old[m_mm] = n_t;
/*
* Calculate the mole numbers of species
*/
// Calculate the mole numbers of species
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
writelogf("START ITERATION %d:\n", iter);
}
/*
* Calculate the mole numbers of species and elements.
*/
// Calculate the mole numbers of species and elements.
double n_t_calc = calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc,
pressureConst);
@ -1086,9 +1029,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
double nCutoff;
bool normalStep = true;
/*
* Decide if we are to do a normal step or a modified step
*/
// Decide if we are to do a normal step or a modified step
size_t iM = npos;
for (m = 0; m < m_mm; m++) {
if (elMoles[m] > 0.001 * elMolesTotal) {
@ -1130,30 +1071,28 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
resid[m_mm] = std::min(resid[m_mm], 0.0);
}
} else {
/*
* Determine whether the matrix should be dumbed down because
* the coefficient matrix of species (with significant concentrations)
* is rank deficient.
*
* The basic idea is that at any time during the calculation only a
* small subset of species with sufficient concentration matters.
* If the rank of the element coefficient matrix for that subset of species
* is less than the number of elements, then the matrix created by
* the Brinkley method below may become singular.
*
* The logic below looks for obvious cases where the current element
* coefficient matrix is rank deficient.
*
* The way around rank-deficiency is to lump-sum the corresponding row
* of the matrix. Note, lump-summing seems to work very well in terms of
* its stability properties, i.e., it heads in the right direction,
* albeit with lousy convergence rates.
*
* NOTE: This probably should be extended to a full blown Gauss-Jordan
* factorization scheme in the future. For Example
* the scheme below would fail for the set: HCl NH4Cl, NH3.
* Hopefully, it's caught by the equal rows logic below.
*/
// Determine whether the matrix should be dumbed down because the
// coefficient matrix of species (with significant concentrations)
// is rank deficient.
//
// The basic idea is that at any time during the calculation only a
// small subset of species with sufficient concentration matters. If
// the rank of the element coefficient matrix for that subset of
// species is less than the number of elements, then the matrix
// created by the Brinkley method below may become singular.
//
// The logic below looks for obvious cases where the current element
// coefficient matrix is rank deficient.
//
// The way around rank-deficiency is to lump-sum the corresponding
// row of the matrix. Note, lump-summing seems to work very well in
// terms of its stability properties, i.e., it heads in the right
// direction, albeit with lousy convergence rates.
//
// NOTE: This probably should be extended to a full blown Gauss-
// Jordan factorization scheme in the future. For Example the scheme
// below would fail for the set: HCl NH4Cl, NH3. Hopefully, it's
// caught by the equal rows logic below.
for (m = 0; m < m_mm; m++) {
lumpSum[m] = 1;
}
@ -1189,9 +1128,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
/*
* Formulate the matrix.
*/
// Formulate the matrix.
for (im = 0; im < m_mm; im++) {
m = m_orderVectorElements[im];
if (im < m_nComponents) {
@ -1212,9 +1149,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
a1(m_mm, m_mm) = 0.0;
/*
* Formulate the residual, resid, and the estimate for the convergence criteria, sum
*/
// Formulate the residual, resid, and the estimate for the
// convergence criteria, sum
sum = 0.0;
for (im = 0; im < m_mm; im++) {
m = m_orderVectorElements[im];
@ -1223,12 +1159,12 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
} else {
resid[m] = 0.0;
}
/*
* For equations with positive and negative coefficients, (electronic charge),
* we must mitigate the convergence criteria by a condition limited by
* finite precision of inverting a matrix.
* Other equations with just positive coefficients aren't limited by this.
*/
// For equations with positive and negative coefficients,
// (electronic charge), we must mitigate the convergence
// criteria by a condition limited by finite precision of
// inverting a matrix. Other equations with just positive
// coefficients aren't limited by this.
if (m == m_eloc) {
tmp = resid[m] / (elMoles[m] + elMolesTotal*1.0E-6 + options.absElemTol);
} else {
@ -1274,21 +1210,17 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
if (DEBUG_MODE_ENABLED && ChemEquil_print_lvl > 0) {
writelogf("(it %d) Convergence = %g\n", iter, sum);
}
/*
* Insist on 20x accuracy compared to the top routine.
* There are instances, for ill-conditioned or
* singular matrices where this is needed to move
* the system to a point where the matrices aren't
* singular.
*/
// Insist on 20x accuracy compared to the top routine. There are
// instances, for ill-conditioned or singular matrices where this is
// needed to move the system to a point where the matrices aren't
// singular.
if (sum < 0.05 * options.relTolerance) {
retn = 0;
break;
}
/*
* Row Sum scaling
*/
// Row Sum scaling
for (m = 0; m <= m_mm; m++) {
tmp = 0.0;
for (n = 0; n <= m_mm; n++) {
@ -1322,27 +1254,26 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
writelogf("] = %10.5g\n", resid[m]);
}
}
/*
* Next Step: We have row-summed the equations.
* However, there are some degenerate cases where two
* rows will be multiplies of each other in terms of
* 0 < m, 0 < m part of the matrix. This occurs on a case
* by case basis, and depends upon the current state of the
* element potential values, which affect the concentrations
* of species.
* So, the way we have found to eliminate this problem is to
* lump-sum one of the rows of the matrix, except for the
* last column, and stick it all on the diagonal.
* Then, we at least have a non-singular matrix, and the
* modified equation moves the corresponding unknown in the
* correct direction.
* The previous row-sum operation has made the identification
* of identical rows much simpler.
*
* Note at least 6E-4 is necessary for the comparison.
* I'm guessing 1.0E-3. If two rows are anywhere close to being
* equivalent, the algorithm can get stuck in an oscillatory mode.
*/
// Next Step: We have row-summed the equations. However, there are
// some degenerate cases where two rows will be multiplies of each
// other in terms of 0 < m, 0 < m part of the matrix. This occurs on
// a case by case basis, and depends upon the current state of the
// element potential values, which affect the concentrations of
// species.
//
// So, the way we have found to eliminate this problem is to lump-
// sum one of the rows of the matrix, except for the last column,
// and stick it all on the diagonal. Then, we at least have a non-
// singular matrix, and the modified equation moves the
// corresponding unknown in the correct direction.
//
// The previous row-sum operation has made the identification of
// identical rows much simpler.
//
// Note at least 6E-4 is necessary for the comparison. I'm guessing
// 1.0E-3. If two rows are anywhere close to being equivalent, the
// algorithm can get stuck in an oscillatory mode.
modifiedMatrix = false;
for (m = 0; m < m_mm; m++) {
size_t sameAsRow = npos;
@ -1401,11 +1332,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
"unused elements.");
}
/*
* Figure out the damping coefficient: Use a delta damping
* coefficient formulation: magnitude of change is capped
* to exp(1).
*/
// Figure out the damping coefficient: Use a delta damping
// coefficient formulation: magnitude of change is capped to exp(1).
beta = 1.0;
for (m = 0; m < m_mm; m++) {
if (resid[m] > 1.0) {
@ -1419,9 +1347,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
writelogf("(it %d) Beta = %g\n", iter, beta);
}
}
/*
* Update the solution vector
*/
// Update the solution vector
for (m = 0; m < m_mm; m++) {
x[m] += beta * resid[m];
}

48
src/equil/MultiPhase.cpp
View file

@ -107,9 +107,8 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
m_np = m_phase.size();
m_nsp += p->nSpecies();
// determine if this phase has new elements
// for each new element, add an entry in the map
// from names to index number + 1:
// determine if this phase has new elements for each new element, add an
// entry in the map from names to index number + 1:
string ename;
// iterate over the elements in this phase
@ -117,10 +116,9 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
for (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_enames.push_back(ename);
@ -135,22 +133,21 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
}
}
// If the mixture temperature hasn't been set, then set the
// temperature and pressure to the values for the phase being
// added. There is no good way to do this. However, this will be overridden later.
// If the mixture temperature hasn't been set, then set the temperature and
// pressure to the values for the phase being added. There is no good way to
// do this. However, this will be overridden later.
if (m_temp == 298.15 && p->temperature() > 2.0E-3) {
m_temp = p->temperature();
m_press = p->pressure();
}
// If this is a solution phase, update the minimum and maximum
// mixture temperatures. Stoichiometric phases are excluded,
// since a mixture may define multiple stoichiometric phases,
// each of which has thermo data valid only over a limited
// range. For example, a mixture might be defined to contain a
// phase representing water ice and one representing liquid
// water, only one of which should be present if the mixture
// represents an equilibrium state.
// If this is a solution phase, update the minimum and maximum mixture
// temperatures. Stoichiometric phases are excluded, since a mixture may
// define multiple stoichiometric phases, each of which has thermo data
// valid only over a limited range. For example, a mixture might be defined
// to contain a phase representing water ice and one representing liquid
// water, only one of which should be present if the mixture represents an
// equilibrium state.
if (p->nSpecies() > 1) {
m_Tmin = std::max(p->minTemp(), m_Tmin);
m_Tmax = std::min(p->maxTemp(), m_Tmax);
@ -172,8 +169,7 @@ void MultiPhase::init()
m_elemAbundances.resize(m_nel, 0.0);
// iterate over the elements
// -> fill in m_atoms(m,k), m_snames(k), m_spphase(k),
// m_sptart(ip)
// -> fill in m_atoms(m,k), m_snames(k), m_spphase(k), m_spstart(ip)
for (m = 0; m < m_nel; m++) {
sym = m_enames[m];
k = 0;
@ -210,8 +206,8 @@ void MultiPhase::init()
}
}
/// set the initial composition within each phase to the
/// mole fractions stored in the phase objects
// set the initial composition within each phase to the
// mole fractions stored in the phase objects
m_init = true;
uploadMoleFractionsFromPhases();
updatePhases();
@ -437,9 +433,7 @@ void MultiPhase::setMolesByName(const std::string& x)
void MultiPhase::getMoles(doublereal* molNum) const
{
/*
* First copy in the mole fractions
*/
// First copy in the mole fractions
copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
size_t ik;
doublereal* dtmp = molNum;
@ -605,8 +599,8 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
Hlow = hnow;
}
} else {
// the current enthalpy is greater than the target; therefore the
// current temperature is too high.
// the current enthalpy is greater than the target;
// therefore the current temperature is too high.
if (m_temp < Thigh) {
Thigh = m_temp;
Hhigh = hnow;

87
src/equil/MultiPhaseEquil.cpp
View file

@ -30,17 +30,15 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
m_incl_element.resize(m_nel_mix,1);
for (m = 0; m < m_nel_mix; m++) {
string enm = mix->elementName(m);
// element 'E' or 'e' represents an electron; this
// requires special handling, so save its index
// for later use
// element 'E' or 'e' represents an electron; this requires special
// handling, so save its index for later use
if (enm == "E" || enm == "e") {
m_eloc = m;
}
// if an element other than electrons is not present in
// the mixture, then exclude it and all species containing
// it from the calculation. Electrons are a special case,
// since a species can have a negative number of 'atoms'
// of electrons (positive ions).
// if an element other than electrons is not present in the mixture,
// then exclude it and all species containing it from the calculation.
// Electrons are a special case, since a species can have a negative
// number of 'atoms' of electrons (positive ions).
if (m_mix->elementMoles(m) <= 0.0 && m != m_eloc) {
m_incl_element[m] = 0;
for (k = 0; k < m_nsp_mix; k++) {
@ -65,17 +63,15 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
}
}
// include pure single-constituent phases only if their thermo
// data are valid for this temperature. This is necessary,
// since some thermo polynomial fits are done only for a
// limited temperature range. For example, using the NASA
// polynomial fits for solid ice and liquid water, if this
// were not done the calculation would predict solid ice to be
// present far above its melting point, since the thermo
// polynomial fits only extend to 273.15 K, and give
// unphysical results above this temperature, leading
// (incorrectly) to Gibbs free energies at high temperature
// lower than for liquid water.
// include pure single-constituent phases only if their thermo data are
// valid for this temperature. This is necessary, since some thermo
// polynomial fits are done only for a limited temperature range. For
// example, using the NASA polynomial fits for solid ice and liquid water,
// if this were not done the calculation would predict solid ice to be
// present far above its melting point, since the thermo polynomial fits
// only extend to 273.15 K, and give unphysical results above this
// temperature, leading (incorrectly) to Gibbs free energies at high
// temperature lower than for liquid water.
size_t ip;
for (k = 0; k < m_nsp_mix; k++) {
ip = m_mix->speciesPhaseIndex(k);
@ -92,8 +88,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
}
}
// Now build the list of all species to be included in the
// calculation.
// Now build the list of all species to be included in the calculation.
for (k = 0; k < m_nsp_mix; k++) {
if (m_incl_species[k] ==1) {
m_nsp++;
@ -131,10 +126,9 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
m_order[k] = k;
}
// if the 'start' flag is set, estimate the initial mole
// numbers by doing a linear Gibbs minimization. In this case,
// only the elemental composition of the initial mixture state
// matters.
// if the 'start' flag is set, estimate the initial mole numbers by doing a
// linear Gibbs minimization. In this case, only the elemental composition
// of the initial mixture state matters.
if (start) {
setInitialMoles(loglevel-1);
}
@ -160,10 +154,9 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
m_force = false;
updateMixMoles();
// At this point, the instance has been created, the species
// to be included have been determined, and an initial
// composition has been selected that has all non-zero mole
// numbers for the included species.
// At this point, the instance has been created, the species to be included
// have been determined, and an initial composition has been selected that
// has all non-zero mole numbers for the included species.
}
doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err,
@ -219,8 +212,7 @@ int MultiPhaseEquil::setInitialMoles(int loglevel)
int iter = 0;
while (redo) {
// choose a set of components based on the current
// composition
// choose a set of components based on the current composition
computeN();
redo = false;
iter++;
@ -331,18 +323,16 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
}
}
// If a pivot is zero, exchange columns. This occurs when
// a species has an elemental composition that is not
// linearly independent of the component species that have
// already been assigned
// If a pivot is zero, exchange columns. This occurs when a species has
// an elemental composition that is not linearly independent of the
// component species that have already been assigned
if (m < nColumns && m_A(m,m) == 0.0) {
// First, we need to find a good candidate for a
// component species to swap in for the one that has
// zero pivot. It must contain element m, be linearly
// independent of the components processed so far
// (m_A(m,k) != 0), and should be a major species if
// possible. We'll choose the species with greatest
// mole fraction that satisfies these criteria.
// First, we need to find a good candidate for a component species
// to swap in for the one that has zero pivot. It must contain
// element m, be linearly independent of the components processed so
// far (m_A(m,k) != 0), and should be a major species if possible.
// We'll choose the species with greatest mole fraction that
// satisfies these criteria.
doublereal maxmoles = -999.0;
size_t kmax = 0;
for (k = m+1; k < nColumns; k++) {
@ -480,10 +470,9 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
} else {
FCTR = 0.9;
}
// if species k is in a multi-species solution phase, then its
// mole number must remain positive, unless the entire phase
// goes away. First we'll determine an upper bound on omega,
// such that all
// if species k is in a multi-species solution phase, then its mole
// number must remain positive, unless the entire phase goes away. First
// we'll determine an upper bound on omega, such that all
if (m_dsoln[k] == 1) {
if ((m_moles[k] > MAJOR_THRESHOLD) || (ik < m_nel)) {
if (m_moles[k] < MAJOR_THRESHOLD) {
@ -516,9 +505,9 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
// now take a step with this scaled omega
step(omegamax, m_work);
// compute the gradient of G at this new position in the
// current direction. If it is positive, then we have overshot
// the minimum. In this case, interpolate back.
// compute the gradient of G at this new position in the current direction.
// If it is positive, then we have overshot the minimum. In this case,
// interpolate back.
doublereal not_mu = 1.0e12;
m_mix->getValidChemPotentials(not_mu, m_mu.data());
doublereal grad1 = 0.0;