Method to solve a pseudo steady state of a nonlinear problem. More...

#include <solveProb.h>

Collaboration diagram for solveProb:

Public Member Functions
	solveProb (ResidEval *resid)
	Constructor for the object.

virtual	~solveProb ()
	Destructor. Deletes the integrator.

int	solve (int ifunc, doublereal time_scale, doublereal reltol)
	Main routine that actually calculates the pseudo steady state of the surface problem.

virtual void	reportState (doublereal *const CSoln) const
	Report the current state of the solution.

virtual void	setBounds (const doublereal botBounds[], const doublereal topBounds[])
	Set the bottom and top bounds on the solution vector.

void	setAtol (const doublereal atol[])

void	setAtolConst (const doublereal atolconst)

Public Attributes
int	m_ioflag

Private Member Functions
	solveProb (const solveProb &right)
	Unimplemented private copy constructor.

solveProb &	operator= (const solveProb &right)
	Unimplemented private assignment operator.

virtual void	print_header (int ioflag, int ifunc, doublereal time_scale, doublereal reltol, doublereal netProdRate[])
	Printing routine that gets called at the start of every invocation.

virtual void	printIteration (int ioflag, doublereal damp, size_t label_d, size_t label_t, doublereal inv_t, doublereal t_real, int iter, doublereal update_norm, doublereal resid_norm, doublereal netProdRate[], doublereal CSolnSP[], doublereal resid[], doublereal wtSpecies[], size_t dim, bool do_time)
	Printing routine that gets called after every iteration.

virtual void	printFinal (int ioflag, doublereal damp, size_t label_d, size_t label_t, doublereal inv_t, doublereal t_real, int iter, doublereal update_norm, doublereal resid_norm, doublereal netProdRateKinSpecies[], const doublereal CSolnSP[], const doublereal resid[], const doublereal wtSpecies[], const doublereal wtRes[], size_t dim, bool do_time)
	Print a summary of the solution.

virtual doublereal	calc_t (doublereal netProdRateSolnSP[], doublereal Csoln[], size_t label, size_t label_old, doublereal *label_factor, int ioflag)
	Calculate a conservative delta T to use in a pseudo-steady state algorithm.

virtual void	calcWeights (doublereal wtSpecies[], doublereal wtResid[], const doublereal CSolnSP[])
	Calculate the solution and residual weights.

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

virtual void	resjac_eval (std::vector< doublereal * > &JacCol, doublereal resid, doublereal CSolnSP, const doublereal *CSolnSPOld, const bool do_time, const doublereal deltaT)
	Main routine that calculates the current residual and Jacobian.

virtual doublereal	calc_damping (doublereal x[], doublereal dxneg[], size_t dim, size_t *label)
	This function calculates a damping factor for the Newton iteration update vector, dxneg, to insure that all solution components stay within prescribed bounds.

Private Attributes
ResidEval *	m_residFunc
	residual function pointer to be solved.

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

vector_fp	m_atol
	m_atol is the absolute tolerance in real units.

doublereal	m_rtol
	m_rtol is the relative error tolerance.

doublereal	m_maxstep
	maximum value of the time step

vector_fp	m_netProductionRatesSave
	Temporary vector with length MAX(1, m_neq)

vector_fp	m_numEqn1
	Temporary vector with length MAX(1, m_neq)

vector_fp	m_numEqn2
	Temporary vector with length MAX(1, m_neq)

vector_fp	m_CSolnSave
	Temporary vector with length MAX(1, m_neq)

vector_fp	m_CSolnSP
	Solution vector.

vector_fp	m_CSolnSPInit
	Saved initial solution vector.

vector_fp	m_CSolnSPOld
	Saved solution vector at the old time step.

vector_fp	m_wtResid
	Weights for the residual norm calculation.

vector_fp	m_wtSpecies
	Weights for the species concentrations norm calculation.

vector_fp	m_resid
	Residual for the surface problem.

vector_int	m_ipiv
	pivots

std::vector< doublereal * >	m_JacCol
	Vector of pointers to the top of the columns of the jacobians.

Array2D	m_Jac
	Jacobian.

vector_fp	m_topBounds
	Top bounds for the solution vector.

vector_fp	m_botBounds
	Bottom bounds for the solution vector.

Detailed Description

Method to solve a pseudo steady state of a nonlinear problem.

The following class handles the solution of a nonlinear problem.

Res_ss(C) = - Res(C) = 0

Optionally a pseudo transient algorithm may be used to relax the residual if it is available.

Res_td(C) = dC/dt - Res(C) = 0;

Res_ss(C) is the steady state residual to be solved. Res_td(C) is the time dependent residual which leads to the steady state residual.

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 the following 4 values:

1: SOLVEPROB_INITIALIZE = This assumes that the initial guess supplied to the routine is far from the correct one. Substantial work plus transient time-stepping is to be expected to find a solution.

2: SOLVEPROB_RESIDUAL = Need to solve the nonlinear problem in order to calculate quantities for a residual calculation (Can expect a moderate change in the solution vector -> try to solve the system by direct methods with no damping first -> then, try time-stepping if the first method fails) A "time_scale" supplied here is used in the algorithm to determine when to shut off time-stepping.

3: SOLVEPROB_JACOBIAN = Calculation of the surface problem is due to the need for a numerical jacobian for the gas-problem. The solution is expected to be very close to the initial guess, and extra accuracy is needed because solution variables have been delta'd from nominal values to create jacobian entries.

4: SOLVEPROB_TRANSIENT = The transient calculation is performed here for an amount of time specified by "time_scale". It is not guaranteed to be time-accurate - just stable and fairly fast. The solution after del_t time is returned, whether it's converged to a steady state or not. This is a poor man's time stepping algorithm.

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 147 of file solveProb.h.

Constructor & Destructor Documentation

solveProb ( ResidEval * resid )

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. Currently, only the default value of BULK_ETCH is supported.

Definition at line 35 of file solveProb.cpp.

References solveProb::m_atol, solveProb::m_botBounds, solveProb::m_CSolnSave, solveProb::m_CSolnSP, solveProb::m_CSolnSPInit, solveProb::m_CSolnSPOld, solveProb::m_ipiv, solveProb::m_Jac, solveProb::m_JacCol, solveProb::m_neq, solveProb::m_netProductionRatesSave, solveProb::m_numEqn1, solveProb::m_numEqn2, solveProb::m_resid, solveProb::m_residFunc, solveProb::m_topBounds, solveProb::m_wtResid, solveProb::m_wtSpecies, ResidEval::nEquations(), Array2D::ptrColumn(), and Array2D::resize().

~solveProb ( )

virtual

Destructor. Deletes the integrator.

Definition at line 72 of file solveProb.cpp.

solveProb ( const solveProb & right )

private

Unimplemented private copy constructor.

Member Function Documentation

solveProb& operator= ( const solveProb & right )

private

Unimplemented private assignment operator.

int solve	(	int	ifunc,
		doublereal	time_scale,
		doublereal	reltol
	)

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.

Parameters

ifunc	Determines the type of solution algorithm to be used. Possible values are SOLVEPROB_INITIALIZE , SOLVEPROB_RESIDUAL SOLVEPROB_JACOBIAN SOLVEPROB_TRANSIENT .
time_scale	Time over which to integrate the surface equations, where applicable
reltol	Relative tolerance to use

Returns: Returns 1 if the surface problem is successfully solved. Returns -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 83 of file solveProb.cpp.

References solveProb::calc_damping(), solveProb::calc_t(), solveProb::calcWeights(), DATA_PTR, solveProb::fun_eval(), ResidEval::getInitialConditions(), solveProb::m_CSolnSP, solveProb::m_CSolnSPInit, solveProb::m_CSolnSPOld, solveProb::m_ipiv, solveProb::m_Jac, solveProb::m_JacCol, solveProb::m_neq, solveProb::m_netProductionRatesSave, solveProb::m_numEqn1, solveProb::m_resid, solveProb::m_residFunc, solveProb::m_wtResid, solveProb::m_wtSpecies, Cantera::npos, solveProb::print_header(), solveProb::printFinal(), solveProb::printIteration(), solveProb::resjac_eval(), clockWC::secondsWC(), and SOLVEPROB_INITIALIZE.

void reportState ( doublereal *const CSoln ) const

virtual

Report the current state of the solution.

Parameters

Report the solution vector for the nonlinear problem

Definition at line 399 of file solveProb.cpp.

References solveProb::m_CSolnSP.

void setBounds	(	const doublereal	botBounds[],
		const doublereal	topBounds[]
	)

virtual

Set the bottom and top bounds on the solution vector.

The default is for the bottom is 0.0, while the default for the top is 1.0

Parameters

botBounds	Vector of bottom bounds
topBounds	vector of top bounds

Definition at line 676 of file solveProb.cpp.

References solveProb::m_botBounds, solveProb::m_neq, and solveProb::m_topBounds.

void print_header	(	int	ioflag,
		int	ifunc,
		doublereal	time_scale,
		doublereal	reltol,
		doublereal	netProdRate[]
	)

privatevirtual

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

Definition at line 700 of file solveProb.cpp.

References InterfaceKinetics::getNetProductionRates(), Phase::id(), solveProb::m_atol, Kinetics::nPhases(), Phase::nSpecies(), SOLVEPROB_INITIALIZE, Phase::speciesName(), Kinetics::surfacePhaseIndex(), and Kinetics::thermo().

Referenced by solveProb::solve().

void printIteration	(	int	ioflag,
		doublereal	damp,
		size_t	label_d,
		size_t	label_t,
		doublereal	inv_t,
		doublereal	t_real,
		int	iter,
		doublereal	update_norm,
		doublereal	resid_norm,
		doublereal	netProdRate[],
		doublereal	CSolnSP[],
		doublereal	resid[],
		doublereal	wtSpecies[],
		size_t	dim,
		bool	do_time
	)

privatevirtual

Printing routine that gets called after every iteration.

Definition at line 797 of file solveProb.cpp.

References DATA_PTR, InterfaceKinetics::getNetProductionRates(), Cantera::int2str(), Kinetics::kineticsSpeciesName(), solveProb::m_numEqn1, and Cantera::npos.

Referenced by solveProb::solve().

void printFinal	(	int	ioflag,
		doublereal	damp,
		size_t	label_d,
		size_t	label_t,
		doublereal	inv_t,
		doublereal	t_real,
		int	iter,
		doublereal	update_norm,
		doublereal	resid_norm,
		doublereal	netProdRateKinSpecies[],
		const doublereal	CSolnSP[],
		const doublereal	resid[],
		const doublereal	wtSpecies[],
		const doublereal	wtRes[],
		size_t	dim,
		bool	do_time
	)

privatevirtual

Print a summary of the solution.

Definition at line 887 of file solveProb.cpp.

References Cantera::int2str(), solveProb::m_neq, solveProb::m_numEqn1, and Cantera::npos.

Referenced by solveProb::solve().

doublereal calc_t	(	doublereal	netProdRateSolnSP[],
		doublereal	Csoln[],
		size_t *	label,
		size_t *	label_old,
		doublereal *	label_factor,
		int	ioflag
	)

privatevirtual

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: Returns the 1. / delta T to be used on the next step

Definition at line 617 of file solveProb.cpp.

References solveProb::m_neq.

Referenced by solveProb::solve().

void calcWeights	(	doublereal	wtSpecies[],
		doublereal	wtResid[],
		const doublereal	CSolnSP[]
	)

privatevirtual

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.
CSolnSP	Solution vector for the surface problem

Definition at line 581 of file solveProb.cpp.

References solveProb::m_atol, solveProb::m_Jac, solveProb::m_neq, and solveProb::m_rtol.

Referenced by solveProb::solve().

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

privatevirtual

Main Function evaluation.

Parameters

resid	output Vector of residuals, length = m_neq
CSolnSP	Vector of species concentrations, unknowns in the problem, length = m_neq
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 414 of file solveProb.cpp.

References solveProb::m_residFunc.

Referenced by solveProb::resjac_eval(), and solveProb::solve().

void resjac_eval	(	std::vector< doublereal * > &	JacCol,
		doublereal *	resid,
		doublereal *	CSolnSP,
		const doublereal *	CSolnSPOld,
		const bool	do_time,
		const doublereal	deltaT
	)

privatevirtual

Main routine that calculates the current residual and Jacobian.

Parameters

JacCol	Vector of pointers to the tops of columns of the 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 432 of file solveProb.cpp.

References DATA_PTR, solveProb::fun_eval(), solveProb::m_atol, solveProb::m_neq, solveProb::m_numEqn2, and ckr::max().

Referenced by solveProb::solve().

doublereal calc_damping	(	doublereal	x[],
		doublereal	dxneg[],
		size_t	dim,
		size_t *	label
	)

privatevirtual

This function calculates a damping factor for the Newton iteration update vector, dxneg, to insure that all solution components stay within prescribed bounds.

The default for this class is that all solution components are bounded between zero and one. this is because the original unknowns were mole fractions and surface site fractions.

dxneg[] = negative of the update vector.

The constant "APPROACH" sets the fraction of the distance to the boundary that the step can take. If the full step would not force any fraction outside of the bounds, then Newton's method is mostly allowed to operate normally. There is also some solution damping employed.

Parameters

x	Vector of the current solution components
dxneg	Vector of the negative of the full solution update vector.
dim	Size of the solution vector
label	return int, stating which solution component caused the most damping.

Definition at line 484 of file solveProb.cpp.

References solveProb::m_atol, solveProb::m_botBounds, solveProb::m_topBounds, ckr::min(), and Cantera::npos.

Referenced by solveProb::solve().

Member Data Documentation

ResidEval* m_residFunc

private

residual function pointer to be solved.

Definition at line 367 of file solveProb.h.

Referenced by solveProb::fun_eval(), solveProb::solve(), and solveProb::solveProb().

size_t m_neq

private

Total number of equations to solve in the implicit problem.

Note, this can be zero, and frequently is

Definition at line 373 of file solveProb.h.

Referenced by solveProb::calc_t(), solveProb::calcWeights(), solveProb::printFinal(), solveProb::resjac_eval(), solveProb::setBounds(), solveProb::solve(), and solveProb::solveProb().

vector_fp m_atol

private

m_atol is the absolute tolerance in real units.

Definition at line 376 of file solveProb.h.

Referenced by solveProb::calc_damping(), solveProb::calcWeights(), solveProb::print_header(), solveProb::resjac_eval(), and solveProb::solveProb().

doublereal m_rtol

private

m_rtol is the relative error tolerance.

Definition at line 379 of file solveProb.h.

Referenced by solveProb::calcWeights().

doublereal m_maxstep

private

maximum value of the time step

units = seconds

Definition at line 385 of file solveProb.h.

vector_fp m_netProductionRatesSave

private

Temporary vector with length MAX(1, m_neq)

Definition at line 388 of file solveProb.h.

Referenced by solveProb::solve(), and solveProb::solveProb().

vector_fp m_numEqn1

private

Temporary vector with length MAX(1, m_neq)

Definition at line 391 of file solveProb.h.

Referenced by solveProb::printFinal(), solveProb::printIteration(), solveProb::solve(), and solveProb::solveProb().

vector_fp m_numEqn2

private

Temporary vector with length MAX(1, m_neq)

Definition at line 394 of file solveProb.h.

Referenced by solveProb::resjac_eval(), and solveProb::solveProb().

vector_fp m_CSolnSave

private

Temporary vector with length MAX(1, m_neq)

Definition at line 397 of file solveProb.h.

Referenced by solveProb::solveProb().

vector_fp m_CSolnSP

private

Solution vector.

length MAX(1, m_neq)

Definition at line 403 of file solveProb.h.

Referenced by solveProb::reportState(), solveProb::solve(), and solveProb::solveProb().

vector_fp m_CSolnSPInit

private

Saved initial solution vector.

length MAX(1, m_neq)

Definition at line 409 of file solveProb.h.

Referenced by solveProb::solve(), and solveProb::solveProb().

vector_fp m_CSolnSPOld

private

Saved solution vector at the old time step.

length MAX(1, m_neq)

Definition at line 415 of file solveProb.h.

Referenced by solveProb::solve(), and solveProb::solveProb().

vector_fp m_wtResid

private

Weights for the residual norm calculation.

length MAX(1, m_neq)

Definition at line 421 of file solveProb.h.

Referenced by solveProb::solve(), and solveProb::solveProb().

vector_fp m_wtSpecies

private

Weights for the species concentrations norm calculation.

length MAX(1, m_neq)

Definition at line 427 of file solveProb.h.

Referenced by solveProb::solve(), and solveProb::solveProb().

vector_fp 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 440 of file solveProb.h.

Referenced by solveProb::solve(), and solveProb::solveProb().

vector_int m_ipiv

private

pivots

length MAX(1, m_neq)

Definition at line 446 of file solveProb.h.

Referenced by solveProb::solve(), and solveProb::solveProb().

std::vector<doublereal*> m_JacCol

private

Vector of pointers to the top of the columns of the jacobians.