Class ChemEquil implements a chemical equilibrium solver for single-phase solutions. More...

#include <ChemEquil.h>

Detailed Description

Class ChemEquil implements a chemical equilibrium solver for single-phase solutions.

It is a "non-stoichiometric" solver in the terminology of Smith and Missen [40], 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 (for example, 2 H2 + O2). In general, if speed is important, this solver should be tried first, and if it fails then use MultiPhaseEquil.

Definition at line 78 of file ChemEquil.h.

Public Member Functions
	ChemEquil (ThermoPhase &s)
	Constructor combined with the initialization function.

int	equilibrate (ThermoPhase &s, const char *XY, int loglevel=0)
	Equilibrate a phase, holding the elemental composition fixed at the initial value found within the ThermoPhase object s.

int	equilibrate (ThermoPhase &s, const char *XY, vector< double > &elMoles, int loglevel=0)
	Compute the equilibrium composition for two specified properties and the specified element moles.

Public Attributes
EquilOpt	options
	Options controlling how the calculation is carried out.

Protected Member Functions
double	nAtoms (size_t k, size_t m) const
	number of atoms of element m in species k.

void	initialize (ThermoPhase &s)
	Prepare for equilibrium calculations.

void	setToEquilState (ThermoPhase &s, const vector< double > &x, double t)
	Set mixture to an equilibrium state consistent with specified element potentials and temperature.

int	setInitialMoles (ThermoPhase &s, vector< double > &elMoleGoal, int loglevel=0)
	Estimate the initial mole numbers.

int	estimateElementPotentials (ThermoPhase &s, vector< double > &lambda, vector< double > &elMolesGoal, int loglevel=0)
	Generate a starting estimate for the element potentials.

int	estimateEP_Brinkley (ThermoPhase &s, vector< double > &lambda, vector< double > &elMoles)
	Do a calculation of the element potentials using the Brinkley method, p.

int	dampStep (ThermoPhase &s, vector< double > &oldx, double oldf, vector< double > &grad, vector< double > &step, vector< double > &x, double &f, vector< double > &elmols, double xval, double yval)
	Find an acceptable step size and take it.

void	equilResidual (ThermoPhase &s, const vector< double > &x, const vector< double > &elmtotal, vector< double > &resid, double xval, double yval, int loglevel=0)
	Evaluates the residual vector F, of length m_mm.

void	equilJacobian (ThermoPhase &s, vector< double > &x, const vector< double > &elmols, DenseMatrix &jac, double xval, double yval, int loglevel=0)

void	adjustEloc (ThermoPhase &s, vector< double > &elMolesGoal)

void	update (const ThermoPhase &s)
	Update internally stored state information.

double	calcEmoles (ThermoPhase &s, vector< double > &x, const double &n_t, const vector< double > &Xmol_i_calc, vector< double > &eMolesCalc, vector< double > &n_i_calc, double pressureConst)
	Given a vector of dimensionless element abundances, this routine calculates the moles of the elements and the moles of the species.

Protected Attributes
ThermoPhase *	m_phase
	Pointer to the ThermoPhase object used to initialize this object.

size_t	m_mm
	number of elements in the phase

size_t	m_kk
	number of species in the phase

size_t	m_skip = npos

size_t	m_nComponents
	This is equal to the rank of the stoichiometric coefficient matrix when it is computed.

function< double(ThermoPhase &)>	m_p1

function< double(ThermoPhase &)>	m_p2

vector< double >	m_molefractions
	Current value of the mole fractions in the single phase. length = m_kk.

double	m_elementTotalSum = 1.0
	Current value of the sum of the element abundances given the current element potentials.

vector< double >	m_elementmolefracs
	Current value of the element mole fractions.

vector< double >	m_reswork

vector< double >	m_jwork1

vector< double >	m_jwork2

vector< double >	m_comp
	Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];.

double	m_temp

double	m_dens

double	m_p0 = OneAtm

size_t	m_eloc = npos
	Index of the element id corresponding to the electric charge of each species.

vector< double >	m_startSoln

vector< double >	m_grt

vector< double >	m_mu_RT

vector< double >	m_muSS_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).

vector< size_t >	m_component

double	m_elemFracCutoff = 1e-100
	element fractional cutoff, below which the element will be zeroed.

bool	m_doResPerturb = false

vector< size_t >	m_orderVectorElements

vector< size_t >	m_orderVectorSpecies

int	m_loglevel
	Verbosity of printed output.

Constructor & Destructor Documentation

◆ ChemEquil()

ChemEquil ( ThermoPhase & s )

Constructor combined with the initialization function.

This constructor initializes the ChemEquil object with everything it needs to start solving equilibrium problems.

Parameters

s	ThermoPhase object that will be used in the equilibrium calls.

Definition at line 43 of file ChemEquil.cpp.

Member Function Documentation

◆ equilibrate() [1/2]

int equilibrate	(	ThermoPhase &	s,
		const char *	XY,
		int	loglevel = `0`
	)

Equilibrate a phase, holding the elemental composition fixed at the initial value found within the ThermoPhase object s.

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.

Definition at line 293 of file ChemEquil.cpp.

◆ equilibrate() [2/2]

int equilibrate	(	ThermoPhase &	s,
		const char *	XY,
		vector< double > &	elMoles,
		int	loglevel = `0`
	)

Compute the equilibrium composition for two specified properties and the specified element moles.

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.

Parameters

s	phase object to be equilibrated
XY	property pair to hold constant
elMoles	specified vector of element abundances.
loglevel	Specify amount of debug logging (0 to disable)

Returns: 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.

Definition at line 301 of file ChemEquil.cpp.

◆ nAtoms()

double nAtoms	(	size_t	k,
		size_t	m
	)		const

inlineprotected

number of atoms of element m in species k.

Definition at line 139 of file ChemEquil.h.

◆ initialize()

void initialize ( ThermoPhase & s )

protected

Prepare for equilibrium calculations.

Parameters

s	object representing the solution phase.

Definition at line 48 of file ChemEquil.cpp.

◆ setToEquilState()

void setToEquilState	(	ThermoPhase &	s,
		const vector< double > &	x,
		double	t
	)

protected

Set mixture to an equilibrium state consistent with specified element potentials and temperature.

Parameters

s	mixture to be updated
x	vector of non-dimensional element potentials \[ \lambda_m/RT \] .
t	temperature in K.

Definition at line 116 of file ChemEquil.cpp.

◆ setInitialMoles()

int setInitialMoles	(	ThermoPhase &	s,
		vector< double > &	elMoleGoal,
		int	loglevel = `0`
	)

protected

Estimate the initial mole numbers.

This version borrows from the MultiPhaseEquil solver.

Definition at line 166 of file ChemEquil.cpp.

◆ estimateElementPotentials()

int estimateElementPotentials	(	ThermoPhase &	s,
		vector< double > &	lambda,
		vector< double > &	elMolesGoal,
		int	loglevel = `0`
	)

protected

Generate a starting estimate for the element potentials.

Definition at line 201 of file ChemEquil.cpp.

◆ estimateEP_Brinkley()

int estimateEP_Brinkley	(	ThermoPhase &	s,
		vector< double > &	lambda,
		vector< double > &	elMoles
	)

protected

Do a calculation of the element potentials using the Brinkley method, p.

129 Smith and Missen [40].

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.

Definition at line 828 of file ChemEquil.cpp.

◆ dampStep()

int dampStep	(	ThermoPhase &	s,
		vector< double > &	oldx,
		double	oldf,
		vector< double > &	grad,
		vector< double > &	step,
		vector< double > &	x,
		double &	f,
		vector< double > &	elmols,
		double	xval,
		double	yval
	)

protected

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.

Definition at line 674 of file ChemEquil.cpp.

◆ equilResidual()

void equilResidual	(	ThermoPhase &	s,
		const vector< double > &	x,
		const vector< double > &	elmtotal,
		vector< double > &	resid,
		double	xval,
		double	yval,
		int	loglevel = `0`
	)

protected

Evaluates the residual vector F, of length m_mm.

Definition at line 713 of file ChemEquil.cpp.

◆ equilJacobian()

void equilJacobian	(	ThermoPhase &	s,
		vector< double > &	x,
		const vector< double > &	elmols,
		DenseMatrix &	jac,
		double	xval,
		double	yval,
		int	loglevel = `0`
	)

protected

Definition at line 759 of file ChemEquil.cpp.

◆ adjustEloc()

void adjustEloc	(	ThermoPhase &	s,
		vector< double > &	elMolesGoal
	)

protected

Definition at line 1297 of file ChemEquil.cpp.

◆ update()

void update ( const ThermoPhase & s )

protected

Update internally stored state information.

Definition at line 137 of file ChemEquil.cpp.

◆ calcEmoles()

double calcEmoles	(	ThermoPhase &	s,
		vector< double > &	x,
		const double &	n_t,
		const vector< double > &	Xmol_i_calc,
		vector< double > &	eMolesCalc,
		vector< double > &	n_i_calc,
		double	pressureConst
	)

protected

Given a vector of dimensionless element abundances, this routine calculates the moles of the elements and the moles of the species.

Parameters

	s	ThermoPhase object
[in]	x	current dimensionless element potentials
[in]	Xmol_i_calc	Mole fractions of the species
[in]	pressureConst	Pressure

Definition at line 792 of file ChemEquil.cpp.

Member Data Documentation

◆ options

EquilOpt options

Options controlling how the calculation is carried out.

See also: EquilOpt

Definition at line 127 of file ChemEquil.h.

◆ m_phase

ThermoPhase* m_phase

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.

Definition at line 136 of file ChemEquil.h.

◆ m_mm

size_t m_mm

protected

number of elements in the phase

Definition at line 251 of file ChemEquil.h.

◆ m_kk

size_t m_kk

protected

number of species in the phase

Definition at line 252 of file ChemEquil.h.

◆ m_skip

size_t m_skip = npos

protected

Definition at line 253 of file ChemEquil.h.

◆ m_nComponents

size_t m_nComponents

protected

This is equal to the rank of the stoichiometric coefficient matrix when it is computed.

It's initialized to m_mm.

Definition at line 257 of file ChemEquil.h.

◆ m_p1

function<double(ThermoPhase&)> m_p1

protected

Definition at line 259 of file ChemEquil.h.

◆ m_p2

function<double(ThermoPhase&)> m_p2

protected

Definition at line 259 of file ChemEquil.h.

◆ m_molefractions

vector<double> m_molefractions

protected

Current value of the mole fractions in the single phase. length = m_kk.

Definition at line 262 of file ChemEquil.h.

◆ m_elementTotalSum

double m_elementTotalSum = 1.0

protected

Current value of the sum of the element abundances given the current element potentials.

Definition at line 266 of file ChemEquil.h.

◆ m_elementmolefracs

vector<double> m_elementmolefracs

protected

Current value of the element mole fractions.

Note these aren't the goal element mole fractions.

Definition at line 270 of file ChemEquil.h.

◆ m_reswork

vector<double> m_reswork

protected

Definition at line 271 of file ChemEquil.h.

◆ m_jwork1

vector<double> m_jwork1

protected

Definition at line 272 of file ChemEquil.h.

◆ m_jwork2

vector<double> m_jwork2

protected

Definition at line 273 of file ChemEquil.h.

◆ m_comp

vector<double> m_comp

protected

Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];.

Definition at line 276 of file ChemEquil.h.

◆ m_temp

double m_temp

protected

Definition at line 277 of file ChemEquil.h.

◆ m_dens

double m_dens

protected

Definition at line 277 of file ChemEquil.h.

◆ m_p0

double m_p0 = OneAtm

protected

Definition at line 278 of file ChemEquil.h.

◆ m_eloc

size_t m_eloc = npos

protected

Index of the element id corresponding to the electric charge of each species.

Equal to -1 if there is no such element id.

Definition at line 282 of file ChemEquil.h.

◆ m_startSoln

vector<double> m_startSoln

protected

Definition at line 284 of file ChemEquil.h.

◆ m_grt

vector<double> m_grt

protected

Definition at line 286 of file ChemEquil.h.

◆ m_mu_RT

vector<double> m_mu_RT

protected

Definition at line 287 of file ChemEquil.h.

◆ m_muSS_RT

vector<double> m_muSS_RT

protected

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).

Definition at line 292 of file ChemEquil.h.

◆ m_component

vector<size_t> m_component

protected

Definition at line 293 of file ChemEquil.h.

◆ m_elemFracCutoff

double m_elemFracCutoff = 1e-100

protected

element fractional cutoff, below which the element will be zeroed.

Definition at line 296 of file ChemEquil.h.

◆ m_doResPerturb

bool m_doResPerturb = false

protected

Definition at line 297 of file ChemEquil.h.

◆ m_orderVectorElements

vector<size_t> m_orderVectorElements

protected

Definition at line 299 of file ChemEquil.h.

◆ m_orderVectorSpecies

vector<size_t> m_orderVectorSpecies

protected

Definition at line 300 of file ChemEquil.h.

◆ m_loglevel

int m_loglevel

protected

Verbosity of printed output.

No messages when m_loglevel == 0. More output as level increases.

Definition at line 304 of file ChemEquil.h.

The documentation for this class was generated from the following files:

Detailed Description

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Constructor & Destructor Documentation

◆ ChemEquil()

Member Function Documentation

◆ equilibrate() [1/2]

◆ equilibrate() [2/2]

◆ nAtoms()

◆ initialize()

◆ setToEquilState()

◆ setInitialMoles()

◆ estimateElementPotentials()

◆ estimateEP_Brinkley()

◆ dampStep()

◆ equilResidual()

◆ equilJacobian()

◆ adjustEloc()

◆ update()

◆ calcEmoles()

Member Data Documentation

◆ options

◆ m_phase

◆ m_mm

◆ m_kk

◆ m_skip

◆ m_nComponents

◆ m_p1

◆ m_p2

◆ m_molefractions

◆ m_elementTotalSum

◆ m_elementmolefracs

◆ m_reswork

◆ m_jwork1

◆ m_jwork2

◆ m_comp

◆ m_temp

◆ m_dens

◆ m_p0

◆ m_eloc

◆ m_startSoln

◆ m_grt

◆ m_mu_RT

◆ m_muSS_RT

◆ m_component

◆ m_elemFracCutoff

◆ m_doResPerturb

◆ m_orderVectorElements

◆ m_orderVectorSpecies

◆ m_loglevel