Cantera  2.1.2
Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
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. More...
 
int solve (int ifunc, doublereal time_scale, doublereal reltol)
 Main routine that actually calculates the pseudo steady state of the surface problem. More...
 
virtual void reportState (doublereal *const CSoln) const
 Report the current state of the solution. More...
 
virtual void setBounds (const doublereal botBounds[], const doublereal topBounds[])
 Set the bottom and top bounds on the solution vector. More...
 
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. More...
 
solveProboperator= (const solveProb &right)
 Unimplemented private assignment operator. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
virtual void calcWeights (doublereal wtSpecies[], doublereal wtResid[], const doublereal CSolnSP[])
 Calculate the solution and residual weights. More...
 
virtual void fun_eval (doublereal *const resid, const doublereal *const CSolnSP, const doublereal *const CSolnSPOld, const bool do_time, const doublereal deltaT)
 Main Function evaluation. More...
 
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. More...
 
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. More...
 

Private Attributes

ResidEvalm_residFunc
 residual function pointer to be solved. More...
 
size_t m_neq
 Total number of equations to solve in the implicit problem. More...
 
vector_fp m_atol
 m_atol is the absolute tolerance in real units. More...
 
doublereal m_rtol
 m_rtol is the relative error tolerance. More...
 
doublereal m_maxstep
 maximum value of the time step More...
 
vector_fp m_netProductionRatesSave
 Temporary vector with length MAX(1, m_neq) More...
 
vector_fp m_numEqn1
 Temporary vector with length MAX(1, m_neq) More...
 
vector_fp m_numEqn2
 Temporary vector with length MAX(1, m_neq) More...
 
vector_fp m_CSolnSave
 Temporary vector with length MAX(1, m_neq) More...
 
vector_fp m_CSolnSP
 Solution vector. More...
 
vector_fp m_CSolnSPInit
 Saved initial solution vector. More...
 
vector_fp m_CSolnSPOld
 Saved solution vector at the old time step. More...
 
vector_fp m_wtResid
 Weights for the residual norm calculation. More...
 
vector_fp m_wtSpecies
 Weights for the species concentrations norm calculation. More...
 
vector_fp m_resid
 Residual for the surface problem. More...
 
vector_int m_ipiv
 pivots More...
 
std::vector< doublereal * > m_JacCol
 Vector of pointers to the top of the columns of the jacobians. More...
 
Array2D m_Jac
 Jacobian. More...
 
vector_fp m_topBounds
 Top bounds for the solution vector. More...
 
vector_fp m_botBounds
 Bottom bounds for the solution vector. More...
 

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

Constructor & Destructor Documentation

solveProb ( ResidEval resid)
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
ifuncDetermines the type of solution algorithm to be used. Possible values are SOLVEPROB_INITIALIZE , SOLVEPROB_RESIDUAL SOLVEPROB_JACOBIAN SOLVEPROB_TRANSIENT .
time_scaleTime over which to integrate the surface equations, where applicable
reltolRelative 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 67 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
[out]CSolnsolution vector for the nonlinear problem

Definition at line 386 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
botBoundsVector of bottom bounds
topBoundsvector of top bounds

Definition at line 609 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
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 723 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 812 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
netProdRateSolnSPOutput variable. Net production rate of all of the species in the solution vector.
Csolnoutput variable. Mole fraction of all of the species in the solution vector
labelOutput variable. Pointer to the value of the species index (kindexSP) that is controlling the time step
label_oldOutput variable. Pointer to the value of the species index (kindexSP) that controlled the time step at the previous iteration
label_factorOutput variable. Pointer to the current factor that is used to indicate the same species is controlling the time step.
ioflagLevel of the output requested.
Returns
Returns the 1. / delta T to be used on the next step

Definition at line 557 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.

Calculate the weighting factors for norms wrt both the species concentration unknowns and the residual unknowns.

Parameters
wtSpeciesWeights to use for the soln unknowns. These are in concentration units
wtResidWeights to sue for the residual unknowns.
CSolnSPSolution vector for the surface problem

Definition at line 532 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  CSolnSPOld,
const bool  do_time,
const doublereal  deltaT 
)
privatevirtual

Main Function evaluation.

This calculates the net production rates of all species

Parameters
residoutput Vector of residuals, length = m_neq
CSolnSPVector of species concentrations, unknowns in the problem, length = m_neq
CSolnSPOldOld Vector of species concentrations, unknowns in the problem, length = m_neq
do_timeCalculate a time dependent residual
deltaTDelta time for time dependent problem.

Definition at line 391 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
JacColVector of pointers to the tops of columns of the Jacobian to be evaluated.
residoutput Vector of residuals, length = m_neq
CSolnSPVector of species concentrations, unknowns in the problem, length = m_neq. These are tweaked in order to derive the columns of the jacobian.
CSolnSPOldOld Vector of species concentrations, unknowns in the problem, length = m_neq
do_timeCalculate a time dependent residual
deltaTDelta time for time dependent problem.

Definition at line 406 of file solveProb.cpp.

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

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
xVector of the current solution components
dxnegVector of the negative of the full solution update vector.
dimSize of the solution vector
labelreturn int, stating which solution component caused the most damping.

Definition at line 440 of file solveProb.cpp.

References solveProb::m_atol, solveProb::m_botBounds, solveProb::m_topBounds, and Cantera::npos.

Referenced by solveProb::solve().

Member Data Documentation

ResidEval* m_residFunc
private

residual function pointer to be solved.

Definition at line 343 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 349 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 352 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 355 of file solveProb.h.

Referenced by solveProb::calcWeights().

doublereal m_maxstep
private

maximum value of the time step

units = seconds

Definition at line 361 of file solveProb.h.

vector_fp m_netProductionRatesSave
private

Temporary vector with length MAX(1, m_neq)

Definition at line 364 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 367 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 370 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 373 of file solveProb.h.

Referenced by solveProb::solveProb().

vector_fp m_CSolnSP
private

Solution vector.

length MAX(1, m_neq)

Definition at line 379 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 385 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 391 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 397 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 403 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 416 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 422 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.

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

Definition at line 429 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 435 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 441 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 447 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: