Method to solve a pseudo steady state surface problem. More...

#include <solveSP.h>

Detailed Description

Method to solve a pseudo steady state surface problem.

The following class handles solving the surface problem. The calculation uses Newton's method to obtain the surface fractions of the surface and bulk species by requiring that the surface species production rate = 0 and that the either the bulk fractions are proportional to their production rates or they are constants.

Currently, the bulk mole fractions are treated as constants. Implementation of their being added to the unknown solution vector is delayed.

Lets introduce the unknown vector for the "surface problem". The surface problem is defined as the evaluation of the surface site fractions for multiple surface phases. The unknown vector will consist of the vector of surface concentrations for each species in each surface vector. Species are grouped first by their surface phases

C_i_j = Concentration of the ith species in the jth surface phase
Nj = number of surface species in the jth surface phase

The unknown solution vector is defined as follows:

C_i_j	kindexSP
C_0_0	0
C_1_0	1
C_2_0	2
. . .	...
C_N0-1_0	N0-1
C_0_1	N0
C_1_1	N0+1
C_2_1	N0+2
. . .	...
C_N1-1_1	NO+N1-1

Note there are a couple of different types of species indices floating around in the formulation of this object.

kindexSP: This is the species index in the contiguous vector of unknowns for the surface problem.

Note, in the future, BULK_DEPOSITION systems will be added, and the solveSP unknown vector will get more complicated. It will include the mole fraction and growth rates of specified bulk phases

Indices which relate to individual kinetics objects use the suffix KSI (kinetics species index).

Solution Method

This routine is typically used within a residual calculation in a large code. It's typically invoked millions of times for large calculations, and it must work every time. Therefore, requirements demand that it be robust but also efficient.

The solution methodology is largely determined by the ifunc parameter, that is input to the solution object. This parameter may have one of the values defined in Surface Problem Solver Methods.

Pseudo time stepping algorithm:

The time step is determined from sdot[], so that the time step doesn't ever change the value of a variable by more than 100%.

This algorithm does use a damped Newton's method to relax the equations. Damping is based on a "delta damping" technique. The solution unknowns are not allowed to vary too much between iterations.

EXTRA_ACCURACY: A constant that is the ratio of the required update norm in this Newton iteration compared to that in the nonlinear solver. A value of 0.1 is used so surface species are safely overconverged.

Functions called:

ct_dgetrf – First half of LAPACK direct solve of a full Matrix
ct_dgetrs – Second half of LAPACK direct solve of a full matrix. Returns solution vector in the right-hand-side vector, resid.

Definition at line 143 of file solveSP.h.

Public Member Functions
	solveSP (ImplicitSurfChem *surfChemPtr, int bulkFunc=BULK_ETCH)
	Constructor for the object. More...

	~solveSP ()=default
	Destructor. More...

	solveSP (const solveSP &)=delete

solveSP &	operator= (const solveSP &)=delete

int	solveSurfProb (int ifunc, double time_scale, double TKelvin, double PGas, double reltol, double abstol)
	Main routine that actually calculates the pseudo steady state of the surface problem. More...

Public Attributes
int	m_ioflag = 0

Private Member Functions
void	print_header (int ioflag, int ifunc, double time_scale, int damping, double reltol, double abstol)
	Printing routine that optionally gets called at the start of every invocation. More...

void	printIteration (int ioflag, double damp, int label_d, int label_t, double inv_t, double t_real, size_t iter, double update_norm, double resid_norm, bool do_time, bool final=false)
	Printing routine that gets called after every iteration. More...

double	calc_t (double netProdRateSolnSP[], double XMolSolnSP[], int label, int label_old, double *label_factor, int ioflag)
	Calculate a conservative delta T to use in a pseudo-steady state algorithm. More...

void	calcWeights (double wtSpecies[], double wtResid[], const Array2D &Jac, const double CSolnSP[], const double abstol, const double reltol)
	Calculate the solution and residual weights. More...

void	updateState (const double *cSurfSpec)
	Update the surface states of the surface phases. More...

void	updateMFSolnSP (double *XMolSolnSP)
	Update mole fraction vector consisting of unknowns in surface problem. More...

void	updateMFKinSpecies (double *XMolKinSp, int isp)
	Update the mole fraction vector for a specific kinetic species vector corresponding to one InterfaceKinetics object. More...

void	evalSurfLarge (const double *CSolnSP)
	Update the vector that keeps track of the largest species in each surface phase. More...

void	fun_eval (double resid, const double CSolnSP, const double *CSolnOldSP, const bool do_time, const double deltaT)
	Main Function evaluation. More...

void	resjac_eval (DenseMatrix &jac, double resid, double CSolnSP, const double *CSolnSPOld, const bool do_time, const double deltaT)
	Main routine that calculates the current residual and Jacobian. More...

Private Attributes
vector< InterfaceKinetics * > &	m_objects
	Vector of interface kinetics objects. More...

size_t	m_neq = 0
	Total number of equations to solve in the implicit problem. More...

int	m_bulkFunc
	This variable determines how the bulk phases are to be handled. More...

size_t	m_numSurfPhases = 0
	Number of surface phases in the surface problem. More...

size_t	m_numTotSurfSpecies = 0
	Total number of surface species in all surface phases. More...

vector< size_t >	m_indexKinObjSurfPhase
	Mapping between the surface phases and the InterfaceKinetics objects. More...

vector< size_t >	m_nSpeciesSurfPhase
	Vector of length number of surface phases containing the number of surface species in each phase. More...

vector< SurfPhase * >	m_ptrsSurfPhase
	Vector of surface phase pointers. More...

vector< size_t >	m_eqnIndexStartSolnPhase
	Index of the start of the unknowns for each solution phase. More...

size_t	m_numBulkPhasesSS = 0
	Total number of volumetric condensed phases included in the steady state problem handled by this routine. More...

vector< size_t >	m_numBulkSpecies
	Vector of number of species in the m_numBulkPhases phases. More...

size_t	m_numTotBulkSpeciesSS = 0
	Total number of species in all bulk phases. More...

vector< ThermoPhase * >	m_bulkPhasePtrs
	Vector of bulk phase pointers, length is equal to m_numBulkPhases. More...

vector< size_t >	m_kinSpecIndex
	Index between the equation index and the position in the kinetic species array for the appropriate kinetics operator. More...

vector< size_t >	m_kinObjIndex
	Index between the equation index and the index of the InterfaceKinetics object. More...

vector< size_t >	m_spSurfLarge
	Vector containing the indices of the largest species in each surface phase. More...

size_t	m_maxTotSpecies = 0
	Maximum number of species in any single kinetics operator -> also maxed wrt the total # of solution species. More...

vector< double >	m_netProductionRatesSave
	Temporary vector with length equal to max m_maxTotSpecies. More...

vector< double >	m_numEqn1
	Temporary vector with length equal to max m_maxTotSpecies. More...

vector< double >	m_numEqn2
	Temporary vector with length equal to max m_maxTotSpecies. More...

vector< double >	m_CSolnSave
	Temporary vector with length equal to max m_maxTotSpecies. More...

vector< double >	m_CSolnSP
	Solution vector. length MAX(1, m_neq) More...

vector< double >	m_CSolnSPInit
	Saved initial solution vector. length MAX(1, m_neq) More...

vector< double >	m_CSolnSPOld
	Saved solution vector at the old time step. length MAX(1, m_neq) More...

vector< double >	m_wtResid
	Weights for the residual norm calculation. length MAX(1, m_neq) More...

vector< double >	m_wtSpecies
	Weights for the species concentrations norm calculation. More...

vector< double >	m_resid
	Residual for the surface problem. More...

vector< double >	m_XMolKinSpecies
	Vector of mole fractions. length m_maxTotSpecies. More...

DenseMatrix	m_Jac
	Jacobian. More...

Constructor & Destructor Documentation

◆ solveSP()

solveSP	(	ImplicitSurfChem *	surfChemPtr,
		int	bulkFunc = `BULK_ETCH`
	)

Constructor for the object.

Parameters

surfChemPtr	Pointer to the ImplicitSurfChem object that defines the surface problem to be solved.
bulkFunc	Integer representing how the bulk phases should be handled. See Surface Problem Bulk Phase Mode. Currently, only the default value of BULK_ETCH is supported.

Definition at line 22 of file solveSP.cpp.

◆ ~solveSP()

~solveSP ( )

default

Destructor.

Member Function Documentation

◆ solveSurfProb()

int solveSurfProb	(	int	ifunc,
		double	time_scale,
		double	TKelvin,
		double	PGas,
		double	reltol,
		double	abstol
	)

Main routine that actually calculates the pseudo steady state of the surface problem.

The actual converged solution is returned as part of the internal state of the InterfaceKinetics objects.

Uses Newton's method to get the surface fractions of the surface and bulk species by requiring that the surface species production rate = 0 and that the bulk fractions are proportional to their production rates.

Parameters

ifunc	Determines the type of solution algorithm to be used. See Surface Problem Solver Methods for possible values.
time_scale	Time over which to integrate the surface equations, where applicable
TKelvin	Temperature (kelvin)
PGas	Pressure (pascals)
reltol	Relative tolerance to use
abstol	absolute tolerance.

Returns: 1 if the surface problem is successfully solved. -1 if the surface problem wasn't solved successfully. Note the actual converged solution is returned as part of the internal state of the InterfaceKinetics objects.

Definition at line 89 of file solveSP.cpp.

◆ print_header()

void print_header	(	int	ioflag,
		int	ifunc,
		double	time_scale,
		int	damping,
		double	reltol,
		double	abstol
	)

private

Printing routine that optionally gets called at the start of every invocation.

Definition at line 602 of file solveSP.cpp.

◆ printIteration()

void printIteration	(	int	ioflag,
		double	damp,
		int	label_d,
		int	label_t,
		double	inv_t,
		double	t_real,
		size_t	iter,
		double	update_norm,
		double	resid_norm,
		bool	do_time,
		bool	final = `false`
	)

private

Printing routine that gets called after every iteration.

Definition at line 651 of file solveSP.cpp.

◆ calc_t()

double calc_t	(	double	netProdRateSolnSP[],
		double	XMolSolnSP[],
		int *	label,
		int *	label_old,
		double *	label_factor,
		int	ioflag
	)

private

Calculate a conservative delta T to use in a pseudo-steady state algorithm.

This routine calculates a pretty conservative 1/del_t based on MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees diagonal dominance.

Small surface fractions are allowed to intervene in the del_t determination, no matter how small. This may be changed. Now minimum changed to 1.0e-12,

Maximum time step set to time_scale.

Parameters

netProdRateSolnSP	Output variable. Net production rate of all of the species in the solution vector.
XMolSolnSP	output variable. Mole fraction of all of the species in the solution vector
label	Output variable. Pointer to the value of the species index (kindexSP) that is controlling the time step
label_old	Output variable. Pointer to the value of the species index (kindexSP) that controlled the time step at the previous iteration
label_factor	Output variable. Pointer to the current factor that is used to indicate the same species is controlling the time step.
ioflag	Level of the output requested.

Returns: the 1. / delta T to be used on the next step

Definition at line 563 of file solveSP.cpp.

◆ calcWeights()

void calcWeights	(	double	wtSpecies[],
		double	wtResid[],
		const Array2D &	Jac,
		const double	CSolnSP[],
		const double	abstol,
		const double	reltol
	)

private

Calculate the solution and residual weights.

Parameters

wtSpecies	Weights to use for the soln unknowns. These are in concentration units
wtResid	Weights to sue for the residual unknowns.
Jac	Jacobian. Row sum scaling is used for the Jacobian
CSolnSP	Solution vector for the surface problem
abstol	Absolute error tolerance
reltol	Relative error tolerance

Definition at line 530 of file solveSP.cpp.

◆ updateState()

void updateState ( const double * cSurfSpec )

private

Update the surface states of the surface phases.

Definition at line 272 of file solveSP.cpp.

◆ updateMFSolnSP()

void updateMFSolnSP ( double * XMolSolnSP )

private

Update mole fraction vector consisting of unknowns in surface problem.

Parameters

XMolSolnSP Vector of mole fractions for the unknowns in the surface problem.

Definition at line 286 of file solveSP.cpp.

◆ updateMFKinSpecies()

void updateMFKinSpecies	(	double *	XMolKinSp,
		int	isp
	)

private

Update the mole fraction vector for a specific kinetic species vector corresponding to one InterfaceKinetics object.

Parameters

XMolKinSp	Mole fraction vector corresponding to a particular kinetic species for a single InterfaceKinetics Object This is a vector over all the species in all of the phases in the InterfaceKinetics object
isp	ID of the InterfaceKinetics Object.

Definition at line 294 of file solveSP.cpp.

◆ evalSurfLarge()

void evalSurfLarge ( const double * CSolnSP )

private

Update the vector that keeps track of the largest species in each surface phase.

Parameters

CSolnSP Vector of the current values of the surface concentrations in all of the surface species.

Definition at line 303 of file solveSP.cpp.

◆ fun_eval()

void fun_eval	(	double *	resid,
		const double *	CSolnSP,
		const double *	CSolnOldSP,
		const bool	do_time,
		const double	deltaT
	)

private

Main Function evaluation.

Parameters

resid	output Vector of residuals, length = m_neq
CSolnSP	Vector of species concentrations, unknowns in the problem, length = m_neq
CSolnOldSP	Old Vector of species concentrations, unknowns in the problem, length = m_neq
do_time	Calculate a time dependent residual
deltaT	Delta time for time dependent problem.

Definition at line 319 of file solveSP.cpp.

◆ resjac_eval()

void resjac_eval	(	DenseMatrix &	jac,
		double *	resid,
		double *	CSolnSP,
		const double *	CSolnSPOld,
		const bool	do_time,
		const double	deltaT
	)

private

Main routine that calculates the current residual and Jacobian.

Parameters

jac	Jacobian to be evaluated.
resid	output Vector of residuals, length = m_neq
CSolnSP	Vector of species concentrations, unknowns in the problem, length = m_neq. These are tweaked in order to derive the columns of the Jacobian.
CSolnSPOld	Old Vector of species concentrations, unknowns in the problem, length = m_neq
do_time	Calculate a time dependent residual
deltaT	Delta time for time dependent problem.

Definition at line 418 of file solveSP.cpp.

Member Data Documentation

◆ m_objects

vector<InterfaceKinetics*>& m_objects

private

Vector of interface kinetics objects.

Each of these is associated with one and only one surface phase.

Definition at line 311 of file solveSP.h.

◆ m_neq

size_t m_neq = 0

private

Total number of equations to solve in the implicit problem.

Note, this can be zero, and frequently is

Definition at line 317 of file solveSP.h.

◆ m_bulkFunc

int m_bulkFunc

private

This variable determines how the bulk phases are to be handled.

Possible values are given in Surface Problem Bulk Phase Mode.

Definition at line 323 of file solveSP.h.

◆ m_numSurfPhases

size_t m_numSurfPhases = 0

private

Number of surface phases in the surface problem.

This number is equal to the number of InterfaceKinetics objects in the problem. (until further noted)

Definition at line 330 of file solveSP.h.

◆ m_numTotSurfSpecies

size_t m_numTotSurfSpecies = 0

private

Total number of surface species in all surface phases.

This is also the number of equations to solve for m_mode=0 system It's equal to the sum of the number of species in each of the m_numSurfPhases.

Definition at line 338 of file solveSP.h.

◆ m_indexKinObjSurfPhase

vector<size_t> m_indexKinObjSurfPhase

private

Mapping between the surface phases and the InterfaceKinetics objects.

Currently this is defined to be a 1-1 mapping (and probably assumed in some places) m_surfKinObjID[i] = i

Definition at line 346 of file solveSP.h.

◆ m_nSpeciesSurfPhase

vector<size_t> m_nSpeciesSurfPhase

private

Vector of length number of surface phases containing the number of surface species in each phase.

Length is equal to the number of surface phases, m_numSurfPhases

Definition at line 353 of file solveSP.h.

◆ m_ptrsSurfPhase

vector<SurfPhase*> m_ptrsSurfPhase

private

Vector of surface phase pointers.

This is created during the constructor Length is equal to the number of surface phases, m_numSurfPhases

Definition at line 360 of file solveSP.h.

◆ m_eqnIndexStartSolnPhase

vector<size_t> m_eqnIndexStartSolnPhase

private

Index of the start of the unknowns for each solution phase.

  i_eqn = m_eqnIndexStartPhase[isp]

isp is the phase id in the list of phases solved by the surface problem.

i_eqn is the equation number of the first unknown in the solution vector corresponding to isp'th phase.

Definition at line 372 of file solveSP.h.

◆ m_numBulkPhasesSS

size_t m_numBulkPhasesSS = 0

private

Total number of volumetric condensed phases included in the steady state problem handled by this routine.

This is equal to or less than the total number of volumetric phases in all of the InterfaceKinetics objects. We usually do not include bulk phases. Bulk phases are only included in the calculation when their domain isn't included in the underlying continuum model conservation equation system.

This is equal to 0, for the time being

Definition at line 385 of file solveSP.h.

◆ m_numBulkSpecies

vector<size_t> m_numBulkSpecies

private

Vector of number of species in the m_numBulkPhases phases.

Length is number of bulk phases

Definition at line 391 of file solveSP.h.

◆ m_numTotBulkSpeciesSS

size_t m_numTotBulkSpeciesSS = 0

private

Total number of species in all bulk phases.

This is also the number of bulk equations to solve when bulk equation solving is turned on.

Definition at line 398 of file solveSP.h.

◆ m_bulkPhasePtrs

vector<ThermoPhase*> m_bulkPhasePtrs

private

Vector of bulk phase pointers, length is equal to m_numBulkPhases.

Definition at line 401 of file solveSP.h.

◆ m_kinSpecIndex

vector<size_t> m_kinSpecIndex

private

Index between the equation index and the position in the kinetic species array for the appropriate kinetics operator.

Length = m_neq.

ksp = m_kinSpecIndex[ieq] ksp is the kinetic species index for the ieq'th equation.

Definition at line 411 of file solveSP.h.

◆ m_kinObjIndex

vector<size_t> m_kinObjIndex

private

Index between the equation index and the index of the InterfaceKinetics object.

Length m_neq

Definition at line 418 of file solveSP.h.

◆ m_spSurfLarge

vector<size_t> m_spSurfLarge

private

Vector containing the indices of the largest species in each surface phase.

k = m_spSurfLarge[i] where k is the local species index, that is, it varies from 0 to (num species in phase - 1) and i is the surface phase index in the problem. Length is equal to m_numSurfPhases.

Definition at line 427 of file solveSP.h.

◆ m_maxTotSpecies

size_t m_maxTotSpecies = 0

private

Maximum number of species in any single kinetics operator -> also maxed wrt the total # of solution species.

Definition at line 431 of file solveSP.h.

◆ m_netProductionRatesSave

vector<double> m_netProductionRatesSave

private

Temporary vector with length equal to max m_maxTotSpecies.

Definition at line 434 of file solveSP.h.

◆ m_numEqn1

vector<double> m_numEqn1

private

Temporary vector with length equal to max m_maxTotSpecies.

Definition at line 437 of file solveSP.h.

◆ m_numEqn2

vector<double> m_numEqn2

private

Temporary vector with length equal to max m_maxTotSpecies.

Definition at line 440 of file solveSP.h.

◆ m_CSolnSave

vector<double> m_CSolnSave

private

Temporary vector with length equal to max m_maxTotSpecies.

Definition at line 443 of file solveSP.h.

◆ m_CSolnSP

vector<double> m_CSolnSP

private

Solution vector. length MAX(1, m_neq)

Definition at line 446 of file solveSP.h.

◆ m_CSolnSPInit

vector<double> m_CSolnSPInit

private

Saved initial solution vector. length MAX(1, m_neq)

Definition at line 449 of file solveSP.h.

◆ m_CSolnSPOld

vector<double> m_CSolnSPOld

private

Saved solution vector at the old time step. length MAX(1, m_neq)

Definition at line 452 of file solveSP.h.

◆ m_wtResid

vector<double> m_wtResid

private

Weights for the residual norm calculation. length MAX(1, m_neq)

Definition at line 455 of file solveSP.h.

◆ m_wtSpecies

vector<double> m_wtSpecies

private

Weights for the species concentrations norm calculation.

length MAX(1, m_neq)

Definition at line 461 of file solveSP.h.

◆ m_resid

vector<double> m_resid

private

Residual for the surface problem.

The residual vector of length "dim" that, that has the value of "sdot" for surface species. The residuals for the bulk species are a function of the sdots for all species in the bulk phase. The last residual of each phase enforces {Sum(fractions) = 1}. After linear solve (dgetrf_ & dgetrs_), resid holds the update vector.

length MAX(1, m_neq)

Definition at line 473 of file solveSP.h.

◆ m_XMolKinSpecies

vector<double> m_XMolKinSpecies

private

Vector of mole fractions. length m_maxTotSpecies.

Definition at line 476 of file solveSP.h.

◆ m_Jac

DenseMatrix m_Jac

private

Jacobian.

m_neq by m_neq computed Jacobian matrix for the local Newton's method.

Definition at line 480 of file solveSP.h.

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

Detailed Description

Solution Method

Pseudo time stepping algorithm:

Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Constructor & Destructor Documentation

◆ solveSP()

◆ ~solveSP()

Member Function Documentation

◆ solveSurfProb()

◆ print_header()

◆ printIteration()

◆ calc_t()

◆ calcWeights()

◆ updateState()

◆ updateMFSolnSP()

◆ updateMFKinSpecies()

◆ evalSurfLarge()

◆ fun_eval()

◆ resjac_eval()

Member Data Documentation

◆ m_objects

◆ m_neq

◆ m_bulkFunc

◆ m_numSurfPhases

◆ m_numTotSurfSpecies

◆ m_indexKinObjSurfPhase

◆ m_nSpeciesSurfPhase

◆ m_ptrsSurfPhase

◆ m_eqnIndexStartSolnPhase

◆ m_numBulkPhasesSS

◆ m_numBulkSpecies

◆ m_numTotBulkSpeciesSS

◆ m_bulkPhasePtrs

◆ m_kinSpecIndex

◆ m_kinObjIndex

◆ m_spSurfLarge

◆ m_maxTotSpecies

◆ m_netProductionRatesSave

◆ m_numEqn1

◆ m_numEqn2

◆ m_CSolnSave

◆ m_CSolnSP

◆ m_CSolnSPInit

◆ m_CSolnSPOld

◆ m_wtResid

◆ m_wtSpecies

◆ m_resid

◆ m_XMolKinSpecies

◆ m_Jac