Cantera  2.0
solveProb Class Reference

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

#include <solveProb.h>

Collaboration diagram for solveProb:
[legend]

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

int m_ioflag

## Private Member Functions

solveProb (const solveProb &right)
Unimplemented private copy constructor.

solveProboperator= (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

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

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

 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.

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.

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.

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.

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.

 vector_fp m_atol
private

m_atol is the absolute tolerance in real units.

Definition at line 376 of file solveProb.h.

 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.

 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 m_JacCol
private

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

The "dim" by "dim" computed Jacobian matrix for the local Newton's method.

Definition at line 454 of file solveProb.h.

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

 Array2D m_Jac
private

Jacobian.

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

Definition at line 461 of file solveProb.h.

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

 vector_fp m_topBounds
private

Top bounds for the solution vector.

This defaults to 1.0

Definition at line 467 of file solveProb.h.

Referenced by solveProb::calc_damping(), solveProb::setBounds(), and solveProb::solveProb().

 vector_fp m_botBounds
private

Bottom bounds for the solution vector.

This defaults to 0.0

Definition at line 473 of file solveProb.h.

Referenced by solveProb::calc_damping(), solveProb::setBounds(), and solveProb::solveProb().

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