research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Cleaned up Doxygen documentation in other "numerics" classes

Browse source
This commit is contained in:
Ray Speth 2013-05-29 00:11:40 +00:00
parent eadf4e3fdd
commit 190c4bf8ce
11 changed files with 195 additions and 1161 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
include/cantera/numerics/DAE_Solver.h
View file

@ -169,7 +169,7 @@ public:
* derivatives of all differential variables.
* @param y Calculated value of the solution vector after the procedure ends
* @param yp Calculated value of the solution derivative after the procedure
* @param The first value of t at which a soluton will be
* @param tout The first value of t at which a soluton will be
* requested (from IDASolve). (This is needed here to
* determine the direction of integration and rough scale
* in the independent variable t.

116
include/cantera/numerics/NonlinearSolver.h
View file

@ -82,12 +82,10 @@ namespace Cantera
*
* Damping is used extensively when relaxing the system.
*
*
* The basic idea is that we predict a direction that is parameterized by an overall coordinate
* value, beta, from zero to one, This may or may not be the same as the value, damp,
* depending upon whether the direction is straight.
*
*
* TIME STEP TYPE
*
* The code solves a nonlinear problem. Frequently the nonlinear problem is created from time-dependent
@ -113,30 +111,19 @@ namespace Cantera
*
* MATRIX SCALING
*
*
*
*
* @code
*
*
* NonlinearSolver *nls = new NonlinearSolver(&r1);
*
* int solnType = NSOLN_TYPE_STEADY_STATE ;
*
* nls->setDeltaBoundsMagnitudes(deltaBounds);
*
* nls->solve_nonlinear_problem(solnType, y_comm, ydot_comm, CJ, time_curr, jac,
* num_newt_its, num_linear_solves, numBacktracks,
* loglevelInput);
*
* @endcode
*
*
* @ingroup numerics
*/
class NonlinearSolver
{
public:
//! Default constructor
/*!
@ -144,23 +131,13 @@ public:
*/
NonlinearSolver(ResidJacEval* func);
//!Copy Constructor for the %ThermoPhase object.
/*!
* @param right Item to be copied
*/
//!Copy Constructor
NonlinearSolver(const NonlinearSolver& right);
//! Destructor
~NonlinearSolver();
//! Assignment operator
/*!
* This is NOT a virtual function.
*
* @param right Reference to %NonlinearSolver object to be
* copied into the
* current one.
*/
NonlinearSolver& operator=(const NonlinearSolver& right);
//! Create solution weights for convergence criteria
@ -176,14 +153,13 @@ public:
*/
void createSolnWeights(const doublereal* const y);
//! L2 norm of the delta of the solution vector
//! L2 norm of the delta of the solution vector
/*!
* calculate the norm of the solution vector. This will
* involve the column scaling of the matrix
*
* The third argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
* The third argument has a default of false. However, if true, then a
* table of the largest values is printed out to standard output.
*
* @param delta_y Vector to take the norm of
* @param title Optional title to be printed out
@ -201,16 +177,14 @@ public:
* Calculate the norm of the residual vector. This may
* involve using the row sum scaling from the matrix problem.
*
* The second argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
* The third argument has a default of false. However, if true, then a
* table of the largest values is printed out to standard output.
*
* @param resid Vector of the residuals
* @param title Optional title to be printed out
* @param printLargest Number of specific entries to be printed
* @param y Current value of y - only used for printouts
*
*
* @return Returns the L2 norm of the delta
*/
doublereal residErrorNorm(const doublereal* const resid, const char* title = 0, const int printLargest = 0,
@ -218,19 +192,18 @@ public:
//! Compute the current residual
/*!
* The current value of the residual is stored in the internal work array m_resid, which is defined
* as mutable
* The current value of the residual is stored in the internal work array
* m_resid, which is defined as mutable
*
* @param time_curr Value of the time
* @param typeCalc Type of the calculation
* @param y_curr Current value of the solution vector
* @param ydot_curr Current value of the time derivative of the solution vector
* @param evalType Base evalulation type
* Defaults to Base_ResidEval
* @param evalType Base evalulation type. Defaults to Base_ResidEval
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
* 1 Means a successful operation.
* 0 or neg value Means an unsuccessful operation.
*/
int doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal* const y_curr,
const doublereal* const ydot_curr,
@ -238,11 +211,9 @@ public:
//! Compute the undamped Newton step
/*!
*
* Compute the undamped Newton step. The residual function is
* evaluated at the current time, t_n, at the current values of the
* solution vector, m_y_n, and the solution time derivative, m_ydot_n.
* The Jacobian is not recomputed.
* The residual function is evaluated at the current time, t_n, at the
* current values of the solution vector, m_y_n, and the solution time
* derivative, m_ydot_n. The Jacobian is not recomputed.
*
* A factored jacobian is reused, if available. If a factored jacobian
* is not available, then the jacobian is factored. Before factoring,
@ -256,15 +227,15 @@ public:
* @param delta_y return value of the raw change in y
* @param jac Jacobian
*
* @return Returns the result code from lapack. A zero means success. Anything
* else indicates a failure.
* @return Returns the result code from lapack. A zero means success.
* Anything else indicates a failure.
*/
int doNewtonSolve(const doublereal time_curr, const doublereal* const y_curr,
const doublereal* const ydot_curr, doublereal* const delta_y,
GeneralMatrix& jac);
//! Compute the newton step, either by direct newton's or by solving a close problem that is represented
//! by a Hessian (
//! Compute the newton step, either by direct newton's or by solving a
//! close problem that is represented by a Hessian
/*!
* This is algorith A.6.5.1 in Dennis / Schnabel
*
@ -290,7 +261,6 @@ public:
* ---------------
* internal m_resid Stored residual is used as input
*
*
* @return Returns the result code from lapack. A zero means success. Anything
* else indicates a failure.
*/
@ -323,7 +293,6 @@ public:
void setDeltaBoundsMagnitudes(const doublereal* const deltaBoundsMagnitudes);
protected:
//! Readjust the trust region vectors
/*!
* The trust region is made up of the trust region vector calculation and the trustDelta_ value
@ -333,7 +302,6 @@ protected:
* The trust region calculate is based on
*
* || delta_x dot 1/trustDeltaX_ || <= trustDelta_
*
*/
void readjustTrustVector();
@ -359,12 +327,9 @@ protected:
*/
doublereal calcTrustDistance(std::vector<doublereal> const& deltaX) const;
public:
//! Bound the step
/*!
*
* Return the factor by which the undamped Newton step 'step0'
* must be multiplied in order to keep all solution components in
* all domains between their specified lower and upper bounds.
@ -397,7 +362,6 @@ public:
//! Set bounds constraints for all variables in the problem
/*!
*
* @param y_low_bounds Vector of lower bounds
* @param y_high_bounds Vector of high bounds
*/
@ -424,8 +388,6 @@ public:
//! Function called to evaluate the jacobian matrix and the current
//! residual vector at the current time step
/*!
*
*
* @param J Jacobian matrix to be filled in
* @param f Right hand side. This routine returns the current
* value of the rhs (output), so that it does
@ -468,13 +430,10 @@ public:
//! Return the factor by which the undamped Newton step 'step0'
//! must be multiplied in order to keep the update within the bounds of an accurate jacobian.
/*!
*
* The idea behind these is that the Jacobian couldn't possibly be representative, if the
* variable is changed by a lot. (true for nonlinear systems, false for linear systems)
* Maximum increase in variable in any one newton iteration:
* factor of 1.5
* Maximum decrease in variable in any one newton iteration:
* factor of 2
* Maximum increase in variable in any one newton iteration: factor of 1.5
* Maximum decrease in variable in any one newton iteration: factor of 2
*
* @param y Initial value of the solution vector
* @param step0 initial proposed step size
@ -507,7 +466,17 @@ public:
* @param writetitle Write a title line
* @param num_backtracks Number of backtracks taken
*
* @return returns an integer indicating what happened.
* @return 1 Successful step was taken: Next step was less than previous step.
* s1 is calculated
* 2 Successful step: Next step's norm is less than 0.8
* 3 Success: The final residual is less than 1.0
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 4 Success: The final residual is less than the residual
* from the previous step.
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 0 Uncertain Success: s1 is about the same as s0
* NSOLN_RETN_FAIL_DAMPSTEP
* Unsuccessful step. We can not find a damping factor that is suitable.
*/
int dampStep(const doublereal time_curr, const doublereal* const y_n_curr,
const doublereal* const ydot_n_curr, doublereal* const step_1,
@ -517,8 +486,7 @@ public:
//! Find the solution to F(X) = 0 by damped Newton iteration.
/*!
* On
* entry, x0 contains an initial estimate of the solution. On
* On entry, x0 contains an initial estimate of the solution. On
* successful return, x1 contains the converged solution.
*
* SolnType = TRANSIENT -> we will assume we are relaxing a transient
@ -538,7 +506,6 @@ public:
* @param num_backtracks Number of backtracking steps taken
* @param loglevelInput Input log level determines the amount of printing.
*
*
* @return A positive value indicates a successful convergence
* -1 Failed convergence
*/
@ -558,7 +525,7 @@ public:
setPreviousTimeStep(const std::vector<doublereal>& y_nm1, const std::vector<doublereal>& ydot_nm1);
private:
//! Set the column scales
//! Set the column scaling vector at the current time
void calcColumnScales();
public:
@ -585,7 +552,6 @@ public:
*/
void setColumnScaling(bool useColScaling, const double* const scaleFactors = 0);
//! Set the rowscaling that are used for the inversion of the matrix
/*!
* Row scaling is set here. Right now the row scaling is set internally in the code.
@ -644,11 +610,8 @@ public:
*/
void getResidWts(doublereal* const residWts) const;
//! Check to see if the nonlinear problem has converged
/*!
*
* @param dampCode Code from the damping routine
* @param s1 Value of the norm of the step change
*
@ -702,6 +665,9 @@ public:
*
* residWeightNorm[i] = residAtol[i] + residRtol * m_rowWtScales[i] / neq
*
* @param residRtol scalar residual relative tolerance
* @param residAtol vector of residual absolute tolerances
*
* @param residNormHandling Parameter that sets the default handling of the residual norms
* 0 The residual weighting vector is calculated to make sure that the solution
* norms are roughly 1 when the residual norm is roughly 1.
@ -710,7 +676,7 @@ public:
* 2 Use the minimum value of the residual weights calculated by method 1 and 2.
* This is the default if this routine is called and this parameter isn't specified.
*/
void setResidualTols(double residRtol, double* residATol, int residNormHandling = 2);
void setResidualTols(double residRtol, double* residAtol, int residNormHandling = 2);
//! Set the value of the maximum # of newton iterations
/*!
@ -760,7 +726,6 @@ public:
*/
void descentComparison(doublereal time_curr ,doublereal* ydot0, doublereal* ydot1, int& numTrials);
//! Setup the parameters for the double dog leg
/*!
* The calls to the doCauchySolve() and doNewtonSolve() routines are done at the main level. This routine comes
@ -806,14 +771,11 @@ public:
* 3 Factor of the first Newton step distance
*
* @param factor Factor to use in combination with the method
*
*/
void setTrustRegionInitializationMethod(int method, doublereal factor);
//! Damp using the dog leg approach
/*!
*
* @param time_curr INPUT Current value of the time
* @param y_n_curr INPUT Current value of the solution vector
* @param ydot_n_curr INPUT Current value of the derivative of the solution vector
@ -857,7 +819,6 @@ public:
* @param ydot_n_1 OUTPUT Time derivatives of solution at the conditions which are evaluated for success
* @param trustDeltaOld INPUT Value of the trust length at the old conditions
*
*
* @return This function returns a code which indicates whether the step will be accepted or not.
* 3 Step passed with flying colors. Try redoing the calculation with a bigger trust region.
* 2 Step didn't pass deltaF requirement. Decrease the size of the next trust region for a retry and return
@ -896,7 +857,6 @@ public:
//! Set the print level from the nonlinear solver
/*!
*
* 0 -> absolutely nothing is printed for a single time step.
* 1 -> One line summary per solve_nonlinear call
* 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
@ -923,8 +883,6 @@ public:
*/
void setSolverScheme(int doDogLeg, int doAffineSolve);
/*
* -----------------------------------------------------------------------------------------------------------------
* MEMBER DATA

20
include/cantera/numerics/ResidEval.h
View file

@ -78,10 +78,8 @@ public:
return (m_alg[k] == 1);
}
/**
* Evaluate the residual function. Called by the
* integrator.
* Evaluate the residual function. Called by the integrator.
* @param t time. (input)
* @param y solution vector. (input)
* @param ydot rate of change of solution vector. (input)
@ -109,11 +107,17 @@ public:
return eval(t, y, DATA_PTR(ydot), r);
}
/**
* Fill the solution and derivative vectors with the initial
* conditions at initial time t0.
* @return 1 Everything is fine
* 0 or neg Something went wrong
//! Fill in the initial conditions
/*!
* Values for both the solution and the value of ydot may be provided.
*
* @param[in] t0 Time
* @param[out] y Solution vector
* @param[out] ydot Rate of change of solution vector.
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
virtual int getInitialConditions(const doublereal t0, doublereal* const y,
doublereal* const ydot) {

51
include/cantera/numerics/ResidJacEval.h
View file

@ -47,41 +47,27 @@ enum ResidEval_Type_Enum {
* A class for full (non-sparse dense matrices with Fortran-compatible data storage.
* The class adds support for identifying what types of calls are made to the residual
* evaluator by adding the ResidEval_Type_Enum class.
*
*
*/
class ResidJacEval : public ResidEval
{
public:
//!Default constructor
/*!
* @param atol Initial value of the global tolerance (defaults to 1.0E-13)
*/
ResidJacEval(doublereal atol = 1.0e-13);
//!Copy Constructor for the %ResidJacEval object
/*!
* @param right Item to be copied
*/
//!Copy Constructor
ResidJacEval(const ResidJacEval& right);
//! Assignment operator
/*!
* This is NOT a virtual function.
*
* @param right Reference to %ResidJacEval object to be copied into the
* current one.
*/
ResidJacEval& operator=(const ResidJacEval& right);
//! Duplication routine for objects which inherit from
//! residJacEval
//! Duplication routine for objects derived from residJacEval
/*!
* This virtual routine can be used to duplicate %ResidJacEval objects
* inherited from %ResidJacEval even if the application only has
* a pointer to %ResidJacEval to work with.
* This virtual routine can be used to duplicate objects which inherit
* from ResidJacEval even if the application only has a pointer to
* ResidJacEval to work with.
*
* These routines are basically wrappers around the derived copy
* constructor.
@ -116,33 +102,10 @@ public:
const int id_x = -1,
const doublereal delta_x = 0.0);
/**
* Evaluate the residual function. Called by the
* integrator.
* @param t time. (input)
* @param y solution vector. (input)
* @param ydot rate of change of solution vector. (input)
* @param r residual vector (output)
*/
virtual int eval(const doublereal t, const doublereal* const y,
const doublereal* const ydot,
doublereal* const r);
//! Fill in the initial conditions
/*!
* Values for both the solution and the value of ydot may be provided.
*
* @param t0 Time (input)
* @param y Solution vector (output)
* @param ydot Rate of change of solution vector. (output)
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
virtual int getInitialConditions(const doublereal t0, doublereal* const y, doublereal* const ydot);
//! Filter the solution predictions
@ -251,7 +214,6 @@ public:
//! This function may be used to create output at various points in the execution of an application.
/*!
*
* @param ifunc identity of the call
* 0 Initial call
* 1 Called at the end of every successful time step
@ -324,7 +286,6 @@ public:
const doublereal* const y, const doublereal* const ydot,
GeneralMatrix& J, doublereal* const resid);
//! Calculate an analytical jacobian and the residual at the current time and values.
/*!
* Only called if the jacFormation method is set to analytical
@ -349,7 +310,6 @@ public:
doublereal* const resid);
protected:
//! constant value of atol
doublereal m_atol;
@ -359,4 +319,3 @@ protected:
}
#endif

34
include/cantera/numerics/RootFind.h
View file

@ -43,18 +43,11 @@ namespace Cantera
//! This means that the rootfinder believes the solution is lower than xmin
#define ROOTFIND_SOLNLOWERTHANXMIN -3
//! This means that the rootfinder believes the solution is higher than xmax
/*!
*
*/
#define ROOTFIND_SOLNHIGHERTHANXMAX -4
//@}
//! Root finder for 1D problems
/*!
*
* The root finder solves a single nonlinear equation described below.
*
* \f[
@ -85,13 +78,10 @@ namespace Cantera
*
* Supplying Hints with General Function Behavior Flags
*
*
*
* Stopping Criteria
*
* Specification of the Stopping Criteria
*
*
* Additional constraints
*
* Bounds Criteria For the Routine
@ -141,34 +131,21 @@ namespace Cantera
*/
class RootFind
{
public:
//! Constructor for the object
/*!
*
* @param resid Pointer to the residual function to be used to calculate f(x)
*/
RootFind(ResidEval* resid);
//! Copy constructor
/*!
* @param r object to be copied
*/
RootFind(const RootFind& r);
//! Destructor. Deletes the integrator.
~RootFind();
//! Assignment operator
/*!
* @param right object to be copied
*
* @return returns a reference to the current object
*/
RootFind& operator=(const RootFind& right);
private:
//! Calculate a deltaX from an input value of x
/*!
@ -216,7 +193,6 @@ private:
bool theSame(doublereal x2, doublereal x1, doublereal factor = 1.0) const;
public:
//! Using a line search method, find the root of a 1D function
/*!
* This routine solves the following equation.
@ -333,7 +309,6 @@ public:
void printTable();
public:
//! Pointer to the residual function evaluator
ResidEval* m_residFunc;
@ -356,7 +331,6 @@ public:
doublereal m_maxstep;
protected:
//! Print level
/*!
* 0 No printing of any kind
@ -368,7 +342,6 @@ protected:
int printLvl;
public:
//! Boolean to turn on the possibility of writing a log file.
bool writeLogAllowed_;
@ -413,7 +386,6 @@ protected:
//! Internal variable tracking f(x) of smallest x tried.
doublereal fx_minTried_;
//! Structure containing the iteration history
struct rfTable {
//@{
@ -428,7 +400,6 @@ protected:
double deltaFConverged;
double delX;
std::string reasoning;
void clear() {
@ -453,16 +424,11 @@ protected:
delX(-1.0E300),
reasoning("") {
};
//@}
};
//! Vector of iteration histories
std::vector<struct rfTable> rfHistory_;
};
}
#endif

50
include/cantera/numerics/solveProb.h
View file

@ -55,13 +55,8 @@ const int SOLVEPROB_RESIDUAL = 2;
const int SOLVEPROB_JACOBIAN = 3;
const int SOLVEPROB_TRANSIENT = 4;
namespace Cantera
{
//! Method to solve a pseudo steady state of a nonlinear problem
/*!
* The following class handles the solution of a nonlinear problem.
@ -76,7 +71,6 @@ namespace Cantera
* 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.
@ -84,11 +78,10 @@ namespace Cantera
* work every time. Therefore, requirements demand that it be robust but also
* efficient.
*
* The solution methodology is largely determined by the <TT>ifunc<\TT> parameter,
* The solution methodology is largely determined by the <TT>ifunc</TT> 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
@ -127,7 +120,6 @@ namespace Cantera
* 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.
@ -146,21 +138,11 @@ namespace Cantera
*/
class solveProb
{
public:
//! Constructor for the object
/*!
* @param surfChemPtr Pointer to the ImplicitSurfChem object that
* defines the surface problem to be solved.
*
* @param bulkFunc Integer representing how the bulk phases
* should be handled. Currently, only the
* default value of BULK_ETCH is supported.
*/
solveProb(ResidEval* resid);
//! Destructor. Deletes the integrator.
virtual ~solveProb();
private:
@ -172,7 +154,6 @@ private:
solveProb& operator=(const solveProb& right);
public:
//! Main routine that actually calculates the pseudo steady state
//! of the surface problem
/*!
@ -197,7 +178,7 @@ public:
//! Report the current state of the solution
/*!
* @param Report the solution vector for the nonlinear problem
* @param[out] CSoln solution vector for the nonlinear problem
*/
virtual void reportState(doublereal* const CSoln) const;
@ -210,20 +191,17 @@ public:
*/
virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[]);
void setAtol(const doublereal atol[]);
void setAtolConst(const doublereal atolconst);
private:
//! Printing routine that gets called at the start of every
//! invocation
//! Printing routine that gets called at the start of every invocation
virtual void print_header(int ioflag, int ifunc, doublereal time_scale,
doublereal reltol,
doublereal netProdRate[]);
#ifdef DEBUG_SOLVEPROB
//! Prints out the residual and Jacobian
virtual void printResJac(int ioflag, int neq, const Array2D& Jac,
doublereal resid[], doublereal wtResid[], doublereal norm);
#endif
@ -236,11 +214,7 @@ private:
doublereal resid[],
doublereal wtSpecies[], size_t dim, bool do_time);
//! Print a summary of the solution
/*!
*
*/
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,
@ -264,7 +238,7 @@ private:
*
* @param netProdRateSolnSP Output variable. Net production rate
* of all of the species in the solution vector.
* @param XMolSolnSP output variable.
* @param Csoln output variable.
* Mole fraction of all of the species in the solution vector
* @param label Output variable. Pointer to the value of the
* species index (kindexSP) that is controlling
@ -286,7 +260,9 @@ private:
//! Calculate the solution and residual weights
/*!
* @param wtSpecies Weights to use for the soln unknowns. These
* Calculate the weighting factors for norms wrt both the species
* concentration unknowns and the residual unknowns.
* @param wtSpecies Weights to use for the soln unknowns. These
* are in concentration units
* @param wtResid Weights to sue for the residual unknowns.
*
@ -313,6 +289,7 @@ private:
//! Main Function evaluation
/*!
* This calculates the net production rates of all species
*
* @param resid output Vector of residuals, length = m_neq
* @param CSolnSP Vector of species concentrations, unknowns in the
@ -323,7 +300,7 @@ private:
* @param deltaT Delta time for time dependent problem.
*/
virtual void fun_eval(doublereal* const resid, const doublereal* const CSolnSP,
const doublereal* const CSolnOldSP, const bool do_time, const doublereal deltaT);
const doublereal* const CSolnSPOld, const bool do_time, const doublereal deltaT);
//! Main routine that calculates the current residual and Jacobian
/*!
@ -445,8 +422,7 @@ private:
*/
vector_int m_ipiv;
//! Vector of pointers to the top of the columns of the
//! jacobians
//! 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.
@ -455,8 +431,7 @@ private:
//! Jacobian
/*!
* m_neq by m_neq computed Jacobian matrix for the
* local Newton's method.
* m_neq by m_neq computed Jacobian matrix for the local Newton's method.
*/
Array2D m_Jac;
@ -472,7 +447,6 @@ private:
*/
vector_fp m_botBounds;
public:
int m_ioflag;
};

590
src/numerics/NonlinearSolver.cpp
View file

@ -1,8 +1,5 @@
/**
*
* @file NonlinearSolver.cpp
*
* Damped Newton solver for 0D and 1D problems
* @file NonlinearSolver.cpp Damped Newton solver for 0D and 1D problems
*/
/*
@ -34,16 +31,16 @@ using namespace std;
namespace Cantera
{
//====================================================================================================================
//-----------------------------------------------------------
// Constants
//-----------------------------------------------------------
//! Dampfactor is the factor by which the damping factor is reduced by when a reduction in step length is warranted
//! Dampfactor is the factor by which the damping factor is reduced by when a
//! reduction in step length is warranted
const doublereal DampFactor = 4.0;
//! Number of damping steps that are carried out before the solution is deemed a failure
//! Number of damping steps that are carried out before the solution is deemed
//! a failure
const int NDAMP = 7;
//====================================================================================================================
//! Print a line of a single repeated character string
/*!
* @param str Character string
@ -80,11 +77,6 @@ bool NonlinearSolver::s_doBothSolvesAndCompare(false);
*/
bool NonlinearSolver::s_alwaysAssumeNewtonGood(false);
//====================================================================================================================
// Default constructor
/*
* @param func Residual and jacobian evaluator function object
*/
NonlinearSolver::NonlinearSolver(ResidJacEval* func) :
m_func(func),
solnType_(NSOLN_TYPE_STEADY_STATE),
@ -206,9 +198,8 @@ NonlinearSolver::NonlinearSolver(ResidJacEval* func) :
deltaX_CP_.resize(neq_, 0.0);
Jd_.resize(neq_, 0.0);
deltaX_trust_.resize(neq_, 1.0);
}
//====================================================================================================================
NonlinearSolver::NonlinearSolver(const NonlinearSolver& right) :
m_func(right.m_func),
solnType_(NSOLN_TYPE_STEADY_STATE),
@ -298,13 +289,12 @@ NonlinearSolver::NonlinearSolver(const NonlinearSolver& right) :
*this =operator=(right);
}
//====================================================================================================================
NonlinearSolver::~NonlinearSolver()
{
delete jacCopyPtr_;
delete HessianPtr_;
}
//====================================================================================================================
NonlinearSolver& NonlinearSolver::operator=(const NonlinearSolver& right)
{
if (this == &right) {
@ -405,18 +395,7 @@ NonlinearSolver& NonlinearSolver::operator=(const NonlinearSolver& right)
return *this;
}
//====================================================================================================================
// Create solution weights for convergence criteria
/*
* We create soln weights from the following formula
*
* wt[i] = rtol * abs(y[i]) + atol[i]
*
* The program always assumes that atol is specific
* to the solution component
*
* @param y vector of the current solution values
*/
void NonlinearSolver::createSolnWeights(const doublereal* const y)
{
for (size_t i = 0; i < neq_; i++) {
@ -428,13 +407,7 @@ void NonlinearSolver::createSolnWeights(const doublereal* const y)
#endif
}
}
//====================================================================================================================
// set bounds constraints for all variables in the problem
/*
*
* @param y_low_bounds Vector of lower bounds
* @param y_high_bounds Vector of high bounds
*/
void NonlinearSolver::setBoundsConstraints(const doublereal* const y_low_bounds,
const doublereal* const y_high_bounds)
{
@ -443,38 +416,23 @@ void NonlinearSolver::setBoundsConstraints(const doublereal* const y_low_bounds,
m_y_high_bounds[i] = y_high_bounds[i];
}
}
//====================================================================================================================
void NonlinearSolver::setSolverScheme(int doDogLeg, int doAffineSolve)
{
doDogLeg_ = doDogLeg;
doAffineSolve_ = doAffineSolve;
}
//====================================================================================================================
std::vector<doublereal> & NonlinearSolver::lowBoundsConstraintVector()
{
return m_y_low_bounds;
}
//====================================================================================================================
std::vector<doublereal> & NonlinearSolver::highBoundsConstraintVector()
{
return m_y_high_bounds;
}
//====================================================================================================================
// L2 norm of the delta of the solution vector
/*
* calculate the norm of the solution vector. This will
* involve the column scaling of the matrix
*
* The third argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
*
* @param delta_y Vector to take the norm of
* @param title Optional title to be printed out
* @param printLargest int indicating how many specific lines should be printed out
* @param dampFactor Current value of the damping factor. Defaults to 1.
* only used for printout out a table.
*/
doublereal NonlinearSolver::solnErrorNorm(const doublereal* const delta_y, const char* title, int printLargest,
const doublereal dampFactor) const
{
@ -551,14 +509,7 @@ doublereal NonlinearSolver::solnErrorNorm(const doublereal* const delta_y, const
}
return sum_norm;
}
//====================================================================================================================
/*
* L2 Norm of the residual
*
* The second argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
*/
doublereal NonlinearSolver::residErrorNorm(const doublereal* const resid, const char* title, const int printLargest,
const doublereal* const y) const
{
@ -639,27 +590,7 @@ doublereal NonlinearSolver::residErrorNorm(const doublereal* const resid, const
}
return sum_norm;
}
//====================================================================================================================
// Set the column scaling that are used for the inversion of the matrix
/*
* There are three ways to do this.
*
* The first method is to set the bool useColScaling to true, leaving the scaling factors unset.
* Then, the column scales will be set to the solution error weighting factors. This has the
* effect of ensuring that all delta variables will have the same order of magnitude at convergence
* end.
*
* The second way is the explicitly set the column factors in the second argument of this function call.
*
* The final way to input the scales is to override the ResidJacEval member function call,
*
* calcSolnScales(double time_n, const double *m_y_n_curr, const double *m_y_nm1, double *m_colScales)
*
* Overriding this function call will trump all other ways to specify the column scaling factors.
*
* @param useColScaling Turn this on if you want to use column scaling in the calculations
* @param scaleFactors A vector of doubles that specifies the column factors.
*/
void NonlinearSolver::setColumnScaling(bool useColScaling, const double* const scaleFactors)
{
if (useColScaling) {
@ -678,23 +609,12 @@ void NonlinearSolver::setColumnScaling(bool useColScaling, const double* const s
m_colScaling = 0;
}
}
//====================================================================================================================
// Set the rowscaling that are used for the inversion of the matrix
/*
* Row scaling is set here. Right now the row scaling is set internally in the code.
*
* @param useRowScaling Turn row scaling on or off.
*/
void NonlinearSolver::setRowScaling(bool useRowScaling)
{
m_rowScaling = useRowScaling;
}
//====================================================================================================================
/*
* calcColumnScales():
*
* Set the column scaling vector at the current time
*/
void NonlinearSolver::calcColumnScales()
{
if (m_colScaling == 1) {
@ -710,18 +630,7 @@ void NonlinearSolver::calcColumnScales()
m_func->calcSolnScales(time_n, DATA_PTR(m_y_n_curr), DATA_PTR(m_y_nm1), DATA_PTR(m_colScales));
}
}
//====================================================================================================================
// Compute the current residual
/*
* @param time_curr Value of the time
* @param typeCalc Type of the calculation
* @param y_curr Current value of the solution vector
* @param ydot_curr Current value of the time derivative of the solution vector
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
int NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal* const y_curr,
const doublereal* const ydot_curr, const ResidEval_Type_Enum evalType) const
{
@ -730,14 +639,7 @@ int NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCa
m_resid_scaled = false;
return retn;
}
//====================================================================================================================
// Scale the matrix
/*
* @param jac Jacobian
* @param y_comm Current value of the solution vector
* @param ydot_comm Current value of the time derivative of the solution vector
* @param time_curr current value of the time
*/
void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, doublereal* const ydot_comm,
doublereal time_curr, int num_newt_its)
{
@ -883,14 +785,7 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
}
}
//====================================================================================================================
// Calculate the scaling factor for translating residual norms into solution norms.
/*
* This routine calls computeResidWts() a couple of times in the calculation of m_ScaleSolnNormToResNorm.
* A more sophisticated routine may do more with signs to get a better value. Perhaps, a series of calculations
* with different signs attached may be in order. Then, m_ScaleSolnNormToResNorm would be calculated
* as the minimum of a series of calculations.
*/
void NonlinearSolver::calcSolnToResNormVector()
{
if (! jacCopyPtr_->factored()) {
@ -949,20 +844,7 @@ void NonlinearSolver::calcSolnToResNormVector()
throw CanteraError("NonlinearSolver::calcSolnToResNormVector()" , "Logic error");
}
}
//====================================================================================================================
// Compute the undamped Newton step based on the current jacobian and an input rhs
/*
* Compute the undamped Newton step. The residual function is
* evaluated at the current time, t_n, at the current values of the
* solution vector, m_y_n_curr, and the solution time derivative, m_ydot_n.
* The Jacobian is not recomputed.
*
* A factored jacobian is reused, if available. If a factored jacobian
* is not available, then the jacobian is factored. Before factoring,
* the jacobian is row and column-scaled. Column scaling is not
* recomputed. The row scales are recomputed here, after column
* scaling has been implemented.
*/
int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal* const y_curr,
const doublereal* const ydot_curr, doublereal* const delta_y,
GeneralMatrix& jac)
@ -1045,22 +927,16 @@ int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal*
m_numLocalLinearSolves++;
return info;
}
//====================================================================================================================
// Compute the newton step, either by direct newton's or by solving a close problem that is represented
// by a Hessian (
/*
* This is algorith A.6.5.1 in Dennis / Schnabel
*
* Compute the QR decomposition
*
* Notes on banded Hessian solve:
* The matrix for jT j has a larger band width. Both the top and bottom band widths
* are doubled, going from KU to KU+KL and KL to KU+KL in size. This is not an impossible increase in cost, but
* has to be considered.
*/
int NonlinearSolver::doAffineNewtonSolve(const doublereal* const y_curr, const doublereal* const ydot_curr,
doublereal* const delta_y, GeneralMatrix& jac)
{
/*
* Notes on banded Hessian solve: The matrix for jT j has a larger band
* width. Both the top and bottom band widths are doubled, going from KU
* to KU+KL and KL to KU+KL in size. This is not an impossible increase
* in cost, but has to be considered.
*/
bool newtonGood = true;
// We can default to QR here ( or not )
jac.useFactorAlgorithm(1);
@ -1397,11 +1273,7 @@ int NonlinearSolver::doAffineNewtonSolve(const doublereal* const y_curr, const
return info;
}
//====================================================================================================================
// Do a steepest descent calculation
/*
* This call must be made on the unfactored jacobian!
*/
doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac)
{
doublereal rowFac = 1.0;
@ -1531,7 +1403,7 @@ doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac)
}
return normSoln;
}
//===================================================================================================================
void NonlinearSolver::descentComparison(doublereal time_curr, doublereal* ydot0, doublereal* ydot1, int& numTrials)
{
doublereal ff = 1.0E-5;
@ -1651,21 +1523,10 @@ void NonlinearSolver::descentComparison(doublereal time_curr, doublereal* ydot0
funcDecreaseNewt2 = 0.5 * (residNewt2 - normResid02) / (ff * sNewt);
printf("\t\t %10.3E %13.5E %22.16E\n", ff, funcDecreaseNewt2, residNewt);
}
}
}
}
//====================================================================================================================
// Setup the parameters for the double dog leg
/*
* The calls to the doCauchySolve() and doNewtonSolve() routines are done at the main level. This routine comes
* after those calls. We calculate the point Nuu_ here, the distances of the dog-legs,
* and the norms of the CP and Newton points in terms of the trust vectors.
*/
void NonlinearSolver::setupDoubleDogleg()
{
/*
@ -1720,16 +1581,8 @@ void NonlinearSolver::setupDoubleDogleg()
*/
normTrust_Newton_ = calcTrustDistance(deltaX_Newton_);
normTrust_CP_ = calcTrustDistance(deltaX_CP_);
}
//====================================================================================================================
// Change the global lambda coordinate into the (leg,alpha) coordinate for the double dogleg
/*
* @param lambda Global value of the distance along the double dogleg
* @param alpha relative value along the particular leg
*
* @return Returns the leg number ( 0, 1, or 2).
*/
int NonlinearSolver::lambdaToLeg(const doublereal lambda, doublereal& alpha) const
{
@ -1743,18 +1596,9 @@ int NonlinearSolver::lambdaToLeg(const doublereal lambda, doublereal& alpha) con
alpha = (lambda * dist_Total_ - dist_R0_ - dist_R1_) / dist_R2_;
return 2;
}
//====================================================================================================================
// Calculated the expected residual along the double dogleg curve.
/*
* @param leg 0, 1, or 2 representing the curves of the dogleg
* @param alpha Relative distance along the particular curve.
*
* @return Returns the expected value of the residual at that point according to the quadratic model.
* The residual at the newton point will always be zero.
*/
doublereal NonlinearSolver::expectedResidLeg(int leg, doublereal alpha) const
{
doublereal resD2, res2, resNorm;
doublereal normResid02 = m_normResid_0 * m_normResid_0 * neq_;
@ -1805,16 +1649,7 @@ doublereal NonlinearSolver::expectedResidLeg(int leg, doublereal alpha) const
return resNorm;
}
//====================================================================================================================
// Here we print out the residual at various points along the double dogleg, comparing against the quadratic model
// in a table format
/*
* @param time_curr INPUT current time
* @param ydot0 INPUT Current value of the derivative of the solution vector for non-time dependent
* determinations
* @param legBest OUTPUT leg of the dogleg that gives the lowest residual
* @param alphaBest OUTPUT distance along dogleg for best result.
*/
void NonlinearSolver::residualComparisonLeg(const doublereal time_curr, const doublereal* const ydot0, int& legBest,
doublereal& alphaBest) const
{
@ -1970,18 +1805,13 @@ void NonlinearSolver::residualComparisonLeg(const doublereal time_curr, const do
}
}
//====================================================================================================================
// Calculate the length of the current trust region in terms of the solution error norm
/*
* We carry out a norm of deltaX_trust_ first. Then, we multiply that value
* by trustDelta_
*/
doublereal NonlinearSolver::trustRegionLength() const
{
norm_deltaX_trust_ = solnErrorNorm(DATA_PTR(deltaX_trust_));
return trustDelta_ * norm_deltaX_trust_;
}
//====================================================================================================================
void NonlinearSolver::setDefaultDeltaBoundsMagnitudes()
{
for (size_t i = 0; i < neq_; i++) {
@ -1989,7 +1819,7 @@ void NonlinearSolver::setDefaultDeltaBoundsMagnitudes()
m_deltaStepMinimum[i] = std::max(m_deltaStepMinimum[i], 0.1 * fabs(m_y_n_curr[i]));
}
}
//====================================================================================================================
void NonlinearSolver::adjustUpStepMinimums()
{
for (size_t i = 0; i < neq_; i++) {
@ -2000,37 +1830,18 @@ void NonlinearSolver::adjustUpStepMinimums()
}
}
//====================================================================================================================
void NonlinearSolver::setDeltaBoundsMagnitudes(const doublereal* const deltaStepMinimum)
{
for (size_t i = 0; i < neq_; i++) {
m_deltaStepMinimum[i] = deltaStepMinimum[i];
}
m_manualDeltaStepSet = 1;
}
//====================================================================================================================
/*
*
* Return the factor by which the undamped Newton step 'step0'
* must be multiplied in order to keep the update within the bounds of an accurate jacobian.
*
* The idea behind these is that the Jacobian couldn't possibly be representative, if the
* variable is changed by a lot. (true for nonlinear systems, false for linear systems)
* Maximum increase in variable in any one newton iteration:
* factor of 1.5
* Maximum decrease in variable in any one newton iteration:
* factor of 2
*
* @param y_n_curr Initial value of the solution vector
* @param step_1 initial proposed step size
*
* @return returns the damping factor
*/
double
NonlinearSolver::deltaBoundStep(const doublereal* const y_n_curr, const doublereal* const step_1)
{
size_t i_fbounds = 0;
int ifbd = 0;
int i_fbd = 0;
@ -2125,13 +1936,7 @@ NonlinearSolver::deltaBoundStep(const doublereal* const y_n_curr, const doublere
return f_delta_bounds;
}
//====================================================================================================================
// Readjust the trust region vectors
/*
* The trust region is made up of the trust region vector calculation and the trustDelta_ value
* We periodically recalculate the trustVector_ values so that they renormalize to the
* correct length.
*/
void NonlinearSolver::readjustTrustVector()
{
doublereal trustDeltaOld = trustDelta_;
@ -2194,11 +1999,7 @@ void NonlinearSolver::readjustTrustVector()
trustNormGoal, deltaXSizeOld, trustDeltaOld, norm_deltaX_trust_, trustDelta_);
}
}
//====================================================================================================================
//! Initialize the size of the trust vector.
/*!
* The algorithm we use is to set it equal to the length of the Distance to the Cauchy point.
*/
void NonlinearSolver::initializeTrustRegion()
{
if (trustRegionInitializationMethod_ == 0) {
@ -2242,16 +2043,6 @@ void NonlinearSolver::initializeTrustRegion()
}
}
//====================================================================================================================
// Fill a dogleg solution step vector
/*
* Previously, we have filled up deltaX_Newton_[], deltaX_CP_[], and Nuu_, so that
* this routine is straightforward.
*
* @param leg Leg of the dog leg you are on (0, 1, or 2)
* @param alpha Relative length along the dog length that you are on.
* @param deltaX Vector to be filled up
*/
void NonlinearSolver::fillDogLegStep(int leg, doublereal alpha, std::vector<doublereal> & deltaX) const
{
if (leg == 0) {
@ -2268,16 +2059,7 @@ void NonlinearSolver::fillDogLegStep(int leg, doublereal alpha, std::vector<doub
}
}
}
//====================================================================================================================
// Calculate the trust distance of a step in the solution variables
/*
* The trust distance is defined as the length of the step according to the norm wrt to the trust region.
* We calculate the trust distance by the following method
*
* trustDist = || delta_x dot 1/trustDeltaX_ || / trustDelta_
*
* @param deltaX Current value of deltaX
*/
doublereal NonlinearSolver::calcTrustDistance(std::vector<doublereal> const& deltaX) const
{
doublereal sum = 0.0;
@ -2288,15 +2070,7 @@ doublereal NonlinearSolver::calcTrustDistance(std::vector<doublereal> const& de
}
return sqrt(sum / neq_) / trustDelta_;
}
//====================================================================================================================
// Given a trust distance, this routine calculates the intersection of the this distance with the
// double dogleg curve
/*
* @param trustDelta (INPUT) Value of the trust distance
* @param lambda (OUTPUT) Returns the internal coordinate of the double dogleg
* @param alpha (OUTPUT) Returns the relative distance along the appropriate leg
* @return leg (OUTPUT) Returns the leg ID (0, 1, or 2)
*/
int NonlinearSolver::calcTrustIntersection(doublereal trustDelta, doublereal& lambda, doublereal& alpha) const
{
doublereal dist;
@ -2335,34 +2109,7 @@ int NonlinearSolver::calcTrustIntersection(doublereal trustDelta, doublereal& la
lambda = dist / dist_Total_;
return 1;
}
//====================================================================================================================
/*
*
* boundStep():
*
* Return the factor by which the undamped Newton step 'step0'
* must be multiplied in order to keep all solution components in
* all domains between their specified lower and upper bounds.
* Other bounds may be applied here as well.
*
* Currently the bounds are hard coded into this routine:
*
* Minimum value for all variables: - 0.01 * m_ewt[i]
* Maximum value = none.
*
* Thus, this means that all solution components are expected
* to be numerical greater than zero in the limit of time step
* truncation errors going to zero.
*
* Delta bounds: The idea behind these is that the Jacobian
* couldn't possibly be representative if the
* variable is changed by a lot. (true for
* nonlinear systems, false for linear systems)
* Maximum increase in variable in any one newton iteration:
* factor of 2
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*/
doublereal NonlinearSolver::boundStep(const doublereal* const y, const doublereal* const step0)
{
size_t i_lower = npos;
@ -2412,30 +2159,7 @@ doublereal NonlinearSolver::boundStep(const doublereal* const y, const doublerea
doublereal f_delta_bounds = deltaBoundStep(y, step0);
return std::min(f_bounds, f_delta_bounds);
}
//===================================================================================================================
// Find a damping coefficient through a look-ahead mechanism
/*
*
* On entry, step0 must contain an undamped Newton step to the
* current solution y0. This method attempts to find a damping coefficient
* such that the next undamped step would have a norm smaller than
* that of step0. If successful, the new solution after taking the
* damped step is returned in y1, and the undamped step at y1 is
* returned in step1.
*
*
* @return 1 Successful step was taken: Next step was less than previous step.
* s1 is calculated
* 2 Successful step: Next step's norm is less than 0.8
* 3 Success: The final residual is less than 1.0
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 4 Success: The final residual is less than the residual
* from the previous step.
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 0 Uncertain Success: s1 is about the same as s0
* NSOLN_RETN_FAIL_DAMPSTEP
* Unsuccessful step. We can not find a damping factor that is suitable.
*/
int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal* const y_n_curr,
const doublereal* const ydot_n_curr, doublereal* const step_1,
doublereal* const y_n_1, doublereal* const ydot_n_1, doublereal* const step_2,
@ -2643,31 +2367,7 @@ int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal* cons
}
return NSOLN_RETN_FAIL_DAMPSTEP;
}
//====================================================================================================================
// Damp using the dog leg approach
/*
*
* @param time_curr INPUT Current value of the time
* @param y_n_curr INPUT Current value of the solution vector
* @param ydot_n_curr INPUT Current value of the derivative of the solution vector
* @param step_1 INPUT First trial step for the first iteration
* @param y_n_1 INPUT First trial value of the solution vector
* @param ydot_n_1 INPUT First trial value of the derivative of the solution vector
* @param s1 OUTPUT Norm of the vector step_1
* @param jac INPUT jacobian
* @param numTrials OUTPUT number of trials taken in the current damping step
*
*
* @return 1 Success: Good step was taken. The predicted residual norm is less than one
* 2 Success: Good step: Next step's norm is less than 0.8
* 3 Success: The final residual is less than 1.0
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 4 Success: The final residual is less than the residual from the previous step.
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 0 Unknown Uncertain Success: s1 is about the same as s0
* NSOLN_RETN_FAIL_DAMPSTEP
* Unsuccessful step. Can not find a damping coefficient that is suitable
*/
int NonlinearSolver::dampDogLeg(const doublereal time_curr, const doublereal* y_n_curr,
const doublereal* ydot_n_curr, std::vector<doublereal> & step_1,
doublereal* const y_n_1, doublereal* const ydot_n_1,
@ -2816,34 +2516,7 @@ int NonlinearSolver::dampDogLeg(const doublereal time_curr, const doublereal* y_
}
return NSOLN_RETN_FAIL_DAMPSTEP;
}
//====================================================================================================================
// Decide whether the current step is acceptable and adjust the trust region size
/*
* This is an extension of algorithm 6.4.5 of Dennis and Schnabel.
*
* Here we decide whether to accept the current step
* At the end of the calculation a new estimate of the trust region is calculated
*
* @param time_curr INPUT Current value of the time
* @param leg INPUT Leg of the dogleg that we are on
* @param alpha INPUT Distance down that leg that we are on
* @param y0 INPUT Current value of the solution vector
* @param ydot0 INPUT Current value of the derivative of the solution vector
* @param step0 INPUT Trial step
* @param y1 OUTPUT Solution values at the conditions which are evaluated for success
* @param ydot1 OUTPUT Time derivatives of solution at the conditions which are evalulated for success
* @param trustDeltaOld INPUT Value of the trust length at the old conditions
*
*
* @return This function returns a code which indicates whether the step will be accepted or not.
* 3 Step passed with flying colors. Try redoing the calculation with a bigger trust region.
* 2 Step didn't pass deltaF requirement. Decrease the size of the next trust region for a retry and return
* 0 The step passed.
* -1 The step size is now too small (||d || < 0.1). A really small step isn't decreasing the function.
* This is an error condition.
* -2 Current value of the solution vector caused a residual error in its evaluation.
* Step is a failure, and the step size must be reduced in order to proceed further.
*/
int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal alpha,
const doublereal* const y_n_curr,
const doublereal* const ydot_n_curr, const std::vector<doublereal> & step_1,
@ -3010,21 +2683,7 @@ int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal
}
return retn;
}
//====================================================================================================================
/*
* solve_nonlinear_problem():
*
* Find the solution to F(X) = 0 by damped Newton iteration. On
* entry, x0 contains an initial estimate of the solution. On
* successful return, x1 contains the converged solution.
*
* SolnType = TRANSIENT -> we will assume we are relaxing a transient
* equation system for now. Will make it more general later,
* if an application comes up.
*
* @return A positive value indicates a successful convergence
* -1 Failed convergence
*/
int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_comm, doublereal* const ydot_comm,
doublereal CJ, doublereal time_curr, GeneralMatrix& jac,
int& num_newt_its, int& num_linear_solves,
@ -3585,36 +3244,14 @@ done:
return retnCode;
}
//====================================================================================================================
//! Set the values for the previous time step
/*!
* We set the values for the previous time step here. These are used in the nonlinear
* solve because they affect the calculation of ydot.
*
* @param y_nm1 Value of the solution vector at the previous time step
* @param ydot_nm1 Value of the solution vector derivative at the previous time step
*/
void NonlinearSolver::
setPreviousTimeStep(const std::vector<doublereal>& y_nm1, const std::vector<doublereal>& ydot_nm1)
{
m_y_nm1 = y_nm1;
m_ydot_nm1 = ydot_nm1;
}
//====================================================================================================================
// Print solution norm contribution
/*
* Prints out the most important entries to the update to the solution vector for the current step
*
* @param step_1 Raw update vector for the current nonlinear step
* @param stepNorm_1 Norm of the vector step_1
* @param step_2 Raw update vector for the next solution value based on the old matrix
* @param stepNorm_2 Norm of the vector step_2
* @param title title of the printout
* @param y_n_curr Old value of the solution
* @param y_n_1 New value of the solution after damping corrections
* @param damp Value of the damping factor
* @param num_entries Number of entries to print out
*/
void NonlinearSolver::
print_solnDelta_norm_contrib(const doublereal* const step_1,
const char* const stepNorm_1,
@ -3702,21 +3339,7 @@ static inline doublereal subtractRD(doublereal a, doublereal b)
}
return diff;
}
//====================================================================================================================
/*
*
* Function called by BEuler to evaluate the Jacobian matrix and the
* current residual at the current time step.
* @param N = The size of the equation system
* @param J = Jacobian matrix to be filled in
* @param f = Right hand side. This routine returns the current
* value of the rhs (output), so that it does
* not have to be computed again.
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* 0 Means an unsuccessful operation
*/
int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f,
doublereal time_curr, doublereal CJ,
doublereal* const y, doublereal* const ydot,
@ -3989,16 +3612,7 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f,
return retn;
}
//====================================================================================================================
// Internal function to calculate the time derivative of the solution at the new step
/*
* Previously, the user must have supplied information about the previous time step for this routine to
* work as intended.
*
* @param order of the BDF method
* @param y_curr current value of the solution
* @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr
*/
void NonlinearSolver::
calc_ydot(const int order, const doublereal* const y_curr, doublereal* const ydot_curr) const
{
@ -4025,49 +3639,20 @@ calc_ydot(const int order, const doublereal* const y_curr, doublereal* const ydo
throw CanteraError("calc_ydot()", "Case not covered");
}
}
//====================================================================================================================
// Apply a filtering process to the new step
/*
* @param timeCurrent Current value of the time
* @param y_current current value of the solution
* @param ydot_current Current value of the solution derivative.
*
* @return Returns the norm of the value of the amount filtered
*/
doublereal NonlinearSolver::filterNewStep(const doublereal timeCurrent,
const doublereal* const ybase, doublereal* const step0)
{
return m_func->filterNewStep(timeCurrent, ybase, step0);
}
//====================================================================================================================
// Apply a filtering process to the new solution
/*
* @param timeCurrent Current value of the time
* @param y_current current value of the solution
* @param ydot_current Current value of the solution derivative.
*
* @return Returns the norm of the value of the amount filtered
*/
doublereal NonlinearSolver::filterNewSolution(const doublereal timeCurrent,
doublereal* const y_current, doublereal* const ydot_current)
{
return m_func->filterSolnPrediction(timeCurrent, y_current);
}
//====================================================================================================================
// Compute the Residual Weights
/*
* The residual weights are defined here to be equal to the inverse of the row scaling factors used to
* row scale the matrix, after column scaling is used. They are multiplied by rtol and an atol factor
* is added as well so that if the residual is less than 1, then the calculation is deemed to be converged.
*
* The basic idea is that a change in the solution vector on the order of the convergence tolerance
* multiplied by [RJC] which is of order one after row scaling should give you the relative weight
* of the row. Values of the residual for that row can then be normalized by the value of this weight.
* When the tolerance in delta x is achieved, the tolerance in the residual should also be achieved
* and should be checked.
*/
void
NonlinearSolver::computeResidWts()
void NonlinearSolver::computeResidWts()
{
ResidWtsReevaluated_ = true;
if (checkUserResidualTols_ == 1) {
@ -4098,34 +3683,15 @@ NonlinearSolver::computeResidWts()
}
}
}
//=====================================================================================================================
// return the residual weights
/*
* @param residWts Vector of length neq_
*/
void
NonlinearSolver::getResidWts(doublereal* const residWts) const
void NonlinearSolver::getResidWts(doublereal* const residWts) const
{
for (size_t i = 0; i < neq_; i++) {
residWts[i] = (m_residWts)[i];
}
}
//=====================================================================================================================
// Check to see if the nonlinear problem has converged
/*
*
* @return integer is returned. If positive, then the problem has converged
* 1 Successful step was taken: Next step's norm is less than 1.0.
* The final residual norm is less than 1.0.
* 2 Successful step: Next step's norm is less than 0.8.
* This step's norm is less than 1.0.
* The residual norm can be anything.
* 3 Success: The final residual is less than 1.0E-2
* The predicted deltaSoln is below 1.0E-2.
* 0 Not converged yet
*/
int
NonlinearSolver::convergenceCheck(int dampCode, doublereal s1)
int NonlinearSolver::convergenceCheck(int dampCode, doublereal s1)
{
int retn = 0;
if (m_dampBound < 0.9999) {
@ -4175,13 +3741,7 @@ NonlinearSolver::convergenceCheck(int dampCode, doublereal s1)
}
return retn;
}
//=====================================================================================================================
// Set the absolute tolerances for the solution variables
/*
* Set the absolute tolerances used in the calculation
*
* @param atol Vector of length neq_ that contains the tolerances to be used for the solution variables
*/
void NonlinearSolver::setAtol(const doublereal* const atol)
{
for (size_t i = 0; i < neq_; i++) {
@ -4192,13 +3752,7 @@ void NonlinearSolver::setAtol(const doublereal* const atol)
atolk_[i]= atol[i];
}
}
//=====================================================================================================================
// Set the relative tolerances for the solution variables
/*
* Set the relative tolerances used in the calculation for the solution variables.
*
* @param rtol single double
*/
void NonlinearSolver::setRtol(const doublereal rtol)
{
if (rtol <= 0.0) {
@ -4207,20 +3761,7 @@ void NonlinearSolver::setRtol(const doublereal rtol)
}
rtol_ = rtol;
}
//=====================================================================================================================
// Set the relative and absolute tolerances for the Residual norm comparisons, if used
/*
*
* residWeightNorm[i] = residAtol[i] + residRtol * m_rowWtScales[i] / neq
*
* @param residNormHandling Parameter that sets the default handling of the residual norms
* 0 The residual weighting vector is calculated to make sure that the solution
* norms are roughly 1 when the residual norm is roughly 1.
* This is the default if this routine is not called.
* 1 Use the user residual norm specified by the parameters in this routine
* 2 Use the minimum value of the residual weights calculated by method 0 and 1.
* This is the default if this routine is called and this parameter isn't specified.
*/
void NonlinearSolver::setResidualTols(double residRtol, double* residATol, int residNormHandling)
{
if (residNormHandling < 0 || residNormHandling > 2) {
@ -4241,11 +3782,10 @@ void NonlinearSolver::setResidualTols(double residRtol, double* residATol, int
}
}
}
//=====================================================================================================================
void NonlinearSolver::setPrintLvl(int printLvl)
{
m_print_flag = printLvl;
}
//=====================================================================================================================
}
}

206
src/numerics/ResidJacEval.cpp
View file

@ -20,23 +20,18 @@ using namespace std;
namespace Cantera
{
//====================================================================================================================
ResidJacEval::ResidJacEval(doublereal atol) :
ResidEval(),
m_atol(atol)
{
}
//====================================================================================================================
// Copy Constructor for the %ResidJacEval object
/*
*/
ResidJacEval::ResidJacEval(const ResidJacEval& right) :
ResidEval()
{
*this = operator=(right);
}
//====================================================================================================================
ResidJacEval& ResidJacEval::operator=(const ResidJacEval& right)
{
if (this == &right) {
@ -50,30 +45,17 @@ ResidJacEval& ResidJacEval::operator=(const ResidJacEval& right)
return *this;
}
//====================================================================================================================
// Duplication routine for objects which inherit from %ResidJacEval
/*
* This virtual routine can be used to duplicate %ResidJacEval objects
* inherited from %ResidJacEval even if the application only has
* a pointer to %ResidJacEval to work with.
*
* These routines are basically wrappers around the derived copy
* constructor.
*/
ResidJacEval* ResidJacEval::duplMyselfAsResidJacEval() const
{
return new ResidJacEval(*this);
}
//====================================================================================================================
int ResidJacEval::nEquations() const
{
return neq_;
}
//====================================================================================================================
// Set a global value of the absolute tolerance
/*
* @param atol Value of atol
*/
void ResidJacEval::setAtol(doublereal atol)
{
m_atol = atol;
@ -82,15 +64,7 @@ void ResidJacEval::setAtol(doublereal atol)
"atol must be greater than zero");
}
}
//====================================================================================================================
//! Fill in the initial conditions
/*!
* Values for both the solution and the value of ydot may be provided.
*
* @param t0 Time (input)
* @param y Solution vector (output)
* @param ydot Rate of change of solution vector. (output)
*/
int ResidJacEval::
getInitialConditions(doublereal t0, doublereal* const y, doublereal* const ydot)
{
@ -104,79 +78,27 @@ getInitialConditions(doublereal t0, doublereal* const y, doublereal* const ydot)
}
return 1;
}
//====================================================================================================================
// This function may be used to create output at various points in the execution of an application.
/*
*
* @param ifunc identity of the call
* 0 Initial call
* 1 Called at the end of every successful time step
* -1 Called at the end of every unsuccessful time step
* 2 Called at the end of every call to integrateRJE()
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input)
*/
void ResidJacEval::
user_out2(const int ifunc, const doublereal t, const doublereal deltaT,
const doublereal* y, const doublereal* ydot)
{
}
//====================================================================================================================
// This function may be used to create output at various points in the execution of an application.
/*
* This routine calls user_out2().
*
* @param ifunc identity of the call
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input)
*/
void ResidJacEval::
user_out(const int ifunc, const doublereal t,
const doublereal* y, const doublereal* ydot)
{
user_out2(ifunc, t, 0.0, y, ydot);
}
//====================================================================================================================
//! Evaluate the time tracking equations, if any
/*!
* Evaluate time integrated quantities that are calculated at the
* end of every successful time step. This call is made once at the end of every successful
* time step that advances the time. It's also made once at the start of the time stepping.
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
*/
int ResidJacEval::
evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal* y,
const doublereal* ydot)
{
return 1;
}
//====================================================================================================================
// Return a vector of delta y's for calculation of the numerical Jacobian
/*
* There is a default algorithm provided.
*
* delta_y[i] = atol[i] + 1.0E-6 ysoln[i]
* delta_y[i] = atol[i] + MAX(1.0E-6 ysoln[i] * 0.01 * solnWeights[i])
*
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param delta_y Value of the delta to be used in calculating the numerical jacobian
* @param solnWeights Value of the solution weights that are used in determining convergence (default = 0)
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* 0 Means an unsuccessful operation
*/
int ResidJacEval::
calcDeltaSolnVariables(const doublereal t, const doublereal* const ySoln,
const doublereal* const ySolnDot, doublereal* const deltaYSoln,
@ -193,16 +115,7 @@ calcDeltaSolnVariables(const doublereal t, const doublereal* const ySoln,
}
return 1;
}
//====================================================================================================================
// Returns a vector of column scale factors that can be used to column scale Jacobians.
/*
* Default to yScales[] = 1.0
*
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param y_old Old Solution vector (input, do not modify)
* @param yScales Value of the column scales
*/
void ResidJacEval::
calcSolnScales(const doublereal t, const doublereal* const ysoln, const doublereal* const ysolnOld,
doublereal* const ysolnScales)
@ -215,47 +128,17 @@ calcSolnScales(const doublereal t, const doublereal* const ysoln, const doublere
}
}
}
//====================================================================================================================
// Filter the solution predictions
/*
* Codes might provide a predicted step change. This routine filters the predicted
* solution vector eliminating illegal directions.
*
* @param t Time (input)
* @param y Solution vector (input, output)
* @param step Proposed step in the solution that will be cropped
*/
doublereal ResidJacEval::filterNewStep(doublereal t, const doublereal* const ybase, doublereal* const step)
{
return 0.0;
}
//====================================================================================================================
// Filter the solution predictions
/*
* Codes might provide a predicted solution vector. This routine filters the predicted
* solution vector.
*
* @param t Time (input)
* @param y Solution vector (input, output)
*/
doublereal ResidJacEval::filterSolnPrediction(doublereal t, doublereal* const y)
{
return 0.0;
}
//====================================================================================================================
// Evaluate any stopping criteria other than a final time limit
/*
* If we are to stop the time integration for any reason other than reaching a final time limit, tout,
* provide a test here. This call is made at the end of every successful time step iteration
*
* @return If true, the the time stepping is stopped. If false, then time stepping is stopped if t >= tout
* Defaults to false.
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
*/
bool ResidJacEval::
evalStoppingCritera(const doublereal t,
const doublereal delta_t,
@ -264,42 +147,13 @@ evalStoppingCritera(const doublereal t,
{
return false;
}
//====================================================================================================================
// Multiply the matrix by another matrix that leads to better conditioning
/*
* Provide a left sided matrix that will multiply the current jacobian, after scaling
* and lead to a better conditioned system.
* This routine is called just before the matrix is factored.
*
* Original Problem:
* J delta_x = - Resid
*
* New problem:
* M (J delta_x) = - M Resid
*
* @param matrix Pointer to the current jacobian (if zero, it's already been factored)
* @param nrows offsets for the matrix
* @param rhs residual vector. This also needs to be lhs multiplied by M
*/
int ResidJacEval::
matrixConditioning(doublereal* const matrix, const int nrows, doublereal* const rhs)
{
return 1;
}
//====================================================================================================================
// Evaluate the residual function
/*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param resid Value of the residual that is computed (output)
* @param evalType Type of the residual being computed (defaults to Base_ResidEval)
* @param id_x Index of the variable that is being numerically differenced to find
* the jacobian (defaults to -1, which indicates that no variable is being
* differenced or that the residual doesn't take this issue into account)
* @param delta_x Value of the delta used in the numerical differencing
*/
int ResidJacEval::
evalResidNJ(const doublereal t, const doublereal deltaT, const doublereal* y,
const doublereal* ydot, doublereal* const resid, const ResidEval_Type_Enum evalType,
@ -308,25 +162,14 @@ evalResidNJ(const doublereal t, const doublereal deltaT, const doublereal* y,
throw CanteraError("ResidJacEval::evalResidNJ()", "Not implemented\n");
return 1;
}
//====================================================================================================================
int ResidJacEval::eval(const doublereal t, const doublereal* const y, const doublereal* const ydot,
doublereal* const r)
{
double deltaT = -1.0;
return evalResidNJ(t, deltaT, y, ydot, r);
}
//====================================================================================================================
// Calculate an analytical jacobian and the residual at the current time and values.
/*
* Only called if the jacFormation method is set to analytical
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param J Reference to the SquareMatrix object to be calculated (output)
* @param resid Value of the residual that is computed (output)
*/
int ResidJacEval::
evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
const doublereal* const y,
@ -337,19 +180,7 @@ evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
doublereal* const* jac_colPts = J.colPts();
return evalJacobianDP(t, delta_t, cj, y, ydot, jac_colPts, resid);
}
//====================================================================================================================
// Calculate an analytical jacobian and the residual at the current time and values.
/*
* Only called if the jacFormation method is set to analytical
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param c_j The current value of the coefficient of the time derivative
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param jac_colPts Reference to the SquareMatrix object to be calculated (output)
* @param resid Value of the residual that is computed (output)
*/
int ResidJacEval::
evalJacobianDP(const doublereal t, const doublereal delta_t,
const doublereal c_j,
@ -361,6 +192,5 @@ evalJacobianDP(const doublereal t, const doublereal delta_t,
throw CanteraError("ResidJacEval::evalJacobianDP()", "Not implemented\n");
return 1;
}
//====================================================================================================================
}

144
src/numerics/RootFind.cpp
View file

@ -40,9 +40,6 @@ namespace Cantera
#define DSIGN(x) (( (x) == (0.0) ) ? (0.0) : ( ((x) > 0.0) ? 1.0 : -1.0 ))
#endif
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
#ifdef DEBUG_MODE
//! Print out a form for the current function evaluation
/*!
@ -64,8 +61,6 @@ static void print_funcEval(FILE* fp, doublereal xval, doublereal fval, int its)
}
#endif
//================================================================================================
// Main constructor
RootFind::RootFind(ResidEval* resid) :
m_residFunc(resid),
m_funcTargetValue(0.0),
@ -88,9 +83,8 @@ RootFind::RootFind(ResidEval* resid) :
x_minTried_(1.0E300),
fx_minTried_(0.0)
{
}
//================================================================================================
RootFind::RootFind(const RootFind& r) :
m_residFunc(r.m_residFunc),
m_funcTargetValue(0.0),
@ -115,12 +109,11 @@ RootFind::RootFind(const RootFind& r) :
{
*this = r;
}
//================================================================================================
// Empty destructor
RootFind::~RootFind()
{
}
//====================================================================================================================
RootFind& RootFind::operator=(const RootFind& right)
{
if (this == &right) {
@ -149,14 +142,7 @@ RootFind& RootFind::operator=(const RootFind& right)
return *this;
}
//================================================================================================
// Calculate a deltaX from an input value of x
/*
* This routine ensure that the deltaX will be greater or equal to DeltaXNorm_
* or 1.0E-14 x
*
* @param x1 input value of x
*/
doublereal RootFind::delXNonzero(doublereal x1) const
{
doublereal deltaX = 1.0E-14 * fabs(x1);
@ -166,14 +152,7 @@ doublereal RootFind::delXNonzero(doublereal x1) const
}
return deltaX;
}
//================================================================================================
// Calculate a deltaX from an input value of x
/*
* This routine ensure that the deltaX will be greater or equal to DeltaXNorm_
* or 1.0E-14 x or deltaXConverged_.
*
* @param x1 input value of x
*/
doublereal RootFind::delXMeaningful(doublereal x1) const
{
doublereal del = delXNonzero(x1);
@ -182,16 +161,7 @@ doublereal RootFind::delXMeaningful(doublereal x1) const
}
return del;
}
//================================================================================================
// Calculate a controlled, nonzero delta between two numbers
/*
* The delta is designed to be greater than or equal to delXNonzero(x) defined above
* with the same sign as the original delta. Therefore if you subtract it from either
* of the two original numbers, you get a different number.
*
* @param x1 first number
* @param x2 second number
*/
double RootFind::deltaXControlled(doublereal x2, doublereal x1) const
{
doublereal sgnn = 1.0;
@ -206,18 +176,7 @@ double RootFind::deltaXControlled(doublereal x2, doublereal x1) const
}
return deltaX;
}
//====================================================================================================================
// Function to decide whether two real numbers are the same or not
/*
* A comparison is made between the two numbers to decide whether they
* are close to one another. This is defined as being within factor * delXMeaningful() of each other.
*
* @param x2 First number
* @param x1 second number
* @param factor Multiplicative factor for delta X. defaults to 1
*
* @return Returns a boolean indicating whether the two numbers are the same or not.
*/
bool RootFind::theSame(doublereal x2, doublereal x1, doublereal factor) const
{
doublereal x = fabs(x2) + fabs(x1);
@ -229,26 +188,13 @@ bool RootFind::theSame(doublereal x2, doublereal x1, doublereal factor) const
}
return false;
}
//====================================================================================================================
/*
* The following calculation is a line search method to find the root of a function
*
*
* xbest Returns the x that satisfies the function
* On input, xbest should contain the best estimate
*
* return:
* 0 Found function
*/
int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal& funcTargetValue, doublereal* xbest)
{
/*
* We store the function target and then actually calculate a modified functional
*
* func = eval(x1) - m_funcTargetValue = 0
*
*
*/
m_funcTargetValue = funcTargetValue;
@ -1176,7 +1122,7 @@ done:
return retn;
}
//====================================================================================================================
doublereal RootFind::func(doublereal x)
{
doublereal r;
@ -1198,15 +1144,7 @@ doublereal RootFind::func(doublereal x)
}
return ff;
}
//====================================================================================================================
// Set the tolerance parameters for the rootfinder
/*
* These tolerance parameters are used on the function value to determine convergence
*
*
* @param rtol Relative tolerance. The default is 10^-5
* @param atol absolute tolerance. The default is 10^-11
*/
void RootFind::setTol(doublereal rtolf, doublereal atolf, doublereal rtolx, doublereal atolx)
{
m_atolf = atolf;
@ -1222,38 +1160,12 @@ void RootFind::setTol(doublereal rtolf, doublereal atolf, doublereal rtolx, doub
m_atolx = atolx;
}
}
//====================================================================================================================
// Set the print level from the rootfinder
/*
*
* 0 -> absolutely nothing is printed for a single time step.
* 1 -> One line summary per solve_nonlinear call
* 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
* 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table
* 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed
* 5 -> Algorithm information on the nonlinear iterates are printed out
* 6 -> Additional info on the nonlinear iterates are printed out
* 7 -> Additional info on the linear solve is printed out.
* 8 -> Info on a per iterate of the linear solve is printed out.
*
* @param printLvl integer value
*/
void RootFind::setPrintLvl(int printlvl)
{
printLvl = printlvl;
}
//====================================================================================================================
// Set the function behavior flag
/*
* If this is true, the function is generally an increasing function of x.
* In particular, if the algorithm is seeking a higher value of f, it will look
* in the positive x direction.
*
* This type of function is needed because this algorithm must deal with regions of f(x) where
* f is not changing with x.
*
* @param value boolean value
*/
void RootFind::setFuncIsGenerallyIncreasing(bool value)
{
if (value) {
@ -1261,18 +1173,7 @@ void RootFind::setFuncIsGenerallyIncreasing(bool value)
}
FuncIsGenerallyIncreasing_ = value;
}
//====================================================================================================================
// Set the function behavior flag
/*
* If this is true, the function is generally a decreasing function of x.
* In particular, if the algorithm is seeking a higher value of f, it will look
* in the negative x direction.
*
* This type of function is needed because this algorithm must deal with regions of f(x) where
* f is not changing with x.
*
* @param value boolean value
*/
void RootFind::setFuncIsGenerallyDecreasing(bool value)
{
if (value) {
@ -1280,33 +1181,19 @@ void RootFind::setFuncIsGenerallyDecreasing(bool value)
}
FuncIsGenerallyDecreasing_ = value;
}
//====================================================================================================================
// Set the nominal value of deltaX
/*
* This sets the value of deltaXNorm_
*
* @param deltaXNorm
*/
void RootFind::setDeltaX(doublereal deltaXNorm)
{
DeltaXnorm_ = deltaXNorm;
specifiedDeltaXnorm_ = 1;
}
//====================================================================================================================
// Set the maximum value of deltaX
/*
* This sets the value of deltaXMax_
*
* @param deltaX
*/
void RootFind::setDeltaXMax(doublereal deltaX)
{
DeltaXMax_ = deltaX;
specifiedDeltaXMax_ = 1;
}
//====================================================================================================================
//====================================================================================================================
void RootFind::printTable()
{
printf("\t----------------------------------------------------------------------------------------------------------------------------------------\n");
@ -1322,6 +1209,5 @@ void RootFind::printTable()
}
printf("\t----------------------------------------------------------------------------------------------------------------------------------------\n");
}
//====================================================================================================================
}

18
src/numerics/funcs.cpp
View file

@ -7,6 +7,8 @@
* See file License.txt for licensing information
*/
#include "cantera/numerics/funcs.h"
#include <vector>
#include <algorithm>
@ -16,7 +18,6 @@ using namespace std;
#include "cantera/base/ct_defs.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/base/stringUtils.h"
#include "cantera/numerics/funcs.h"
#include "cantera/numerics/polyfit.h"
#ifndef FTN_TRAILING_UNDERSCORE
@ -38,21 +39,6 @@ extern "C" {
namespace Cantera
{
// Linearly interpolate a function defined on a discrete grid.
/*
* Vector xpts contains a monotonic sequence of grid points, and
* vector fpts contains function values defined at these points.
* The value returned is the linear interpolate at point x.
* If x is outside the range of xpts, the value of fpts at the
* nearest end is returned.
*
* @param x value of the x coordinate
* @param xpts value of the grid points
* @param fpts value of the interpolant at the grid points
*
* @return Returned value is the value of of the interpolated
* function at x.
*/
doublereal linearInterp(doublereal x, const vector_fp& xpts,
const vector_fp& fpts)
{

125
src/numerics/solveProb.cpp
View file

@ -30,8 +30,6 @@ namespace Cantera
static doublereal calcWeightedNorm(const doublereal [], const doublereal dx[], size_t);
//================================================================================================
// Main constructor
solveProb::solveProb(ResidEval* resid) :
m_residFunc(resid),
m_neq(0),
@ -67,22 +65,20 @@ solveProb::solveProb(ResidEval* resid) :
}
}
//================================================================================================
// Empty destructor
solveProb::~solveProb()
{
}
//================================================================================================
/*
* The following calculation is a 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.
*/
int solveProb::solve(int ifunc, doublereal time_scale,
doublereal reltol)
{
/*
* The following calculation is a 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.
*/
doublereal EXTRA_ACCURACY = 0.001;
if (ifunc == SOLVEPROB_JACOBIAN) {
EXTRA_ACCURACY *= 0.001;
@ -392,43 +388,27 @@ int solveProb::solve(int ifunc, doublereal time_scale,
}
return 0;
}
//================================================================================================
/*
* Update the surface states of the surface phases.
*/
void solveProb::reportState(doublereal* const CSolnSP) const
{
std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), CSolnSP);
}
//================================================================================================
/*
* This calculates the net production rates of all species
*
* This calculates the function eval.
* (should switch to special_species formulation for sum condition)
*
* @internal
* This routine uses the m_numEqn1 and m_netProductionRatesSave vectors
* as temporary internal storage.
*/
void solveProb::fun_eval(doublereal* const resid, const doublereal* const CSoln,
const doublereal* const CSolnOld, const bool do_time,
const doublereal deltaT)
{
/*
* This routine uses the m_numEqn1 and m_netProductionRatesSave vectors
* as temporary internal storage.
*/
if (do_time) {
m_residFunc->evalSimpleTD(0.0, CSoln, CSolnOld, deltaT, resid);
} else {
m_residFunc->evalSS(0.0, CSoln, resid);
}
}
//================================================================================================
/*
* Calculate the Jacobian and residual
*
* @internal
* This routine uses the m_numEqn2 vector
* as temporary internal storage.
*/
void solveProb::resjac_eval(std::vector<doublereal*> &JacCol,
doublereal resid[], doublereal CSoln[],
const doublereal CSolnOld[], const bool do_time,
@ -452,6 +432,7 @@ void solveProb::resjac_eval(std::vector<doublereal*> &JacCol,
}
dc = std::max(1.0E-11 * sd, fabs(cSave) * 1.0E-6);
CSoln[kCol] += dc;
// Use the m_numEqn2 vector as temporary internal storage.
fun_eval(DATA_PTR(m_numEqn2), CSoln, CSolnOld, do_time, deltaT);
col_j = JacCol[kCol];
for (size_t i = 0; i < m_neq; i++) {
@ -461,26 +442,9 @@ void solveProb::resjac_eval(std::vector<doublereal*> &JacCol,
}
}
//================================================================================================
#define APPROACH 0.50
// 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.
*
* @param x Vector of the current solution components
* @param dxneg Vector of the negative of the full solution update vector.
* @param dim Size of the solution vector
* @param label return int, stating which solution component caused the most damping.
*/
doublereal solveProb::calc_damping(doublereal x[], doublereal dxneg[], size_t dim, size_t* label)
{
doublereal damp = 1.0, xnew, xtop, xbot;
@ -554,7 +518,7 @@ doublereal solveProb::calc_damping(doublereal x[], doublereal dxneg[], size_t di
}
#undef APPROACH
//================================================================================================
/*
* This function calculates the norm of an update, dx[],
* based on the weighted values of x.
@ -572,12 +536,7 @@ static doublereal calcWeightedNorm(const doublereal wtX[], const doublereal dx[]
}
return sqrt(norm/dim);
}
//================================================================================================
/*
* Calculate the weighting factors for norms wrt both the species
* concentration unknowns and the residual unknowns.
*
*/
void solveProb::calcWeights(doublereal wtSpecies[], doublereal wtResid[],
const doublereal CSoln[])
{
@ -601,18 +560,7 @@ void solveProb::calcWeights(doublereal wtSpecies[], doublereal wtResid[],
}
}
}
//================================================================================================
/*
* 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.
*/
doublereal solveProb::
calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[],
size_t* label, size_t* label_old, doublereal* label_factor, int ioflag)
@ -665,14 +613,7 @@ calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[],
return inv_timeScale;
}
//====================================================================================================================
// 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
*
* @param botBounds Vector of bottom bounds
* @param topBounds vector of top bounds
*/
void solveProb::setBounds(const doublereal botBounds[], const doublereal topBounds[])
{
for (size_t k = 0; k < m_neq; k++) {
@ -680,11 +621,7 @@ void solveProb::setBounds(const doublereal botBounds[], const doublereal topBoun
m_topBounds[k] = topBounds[k];
}
}
//====================================================================================================================
/*
* printResJac(): prints out the residual and Jacobian.
*
*/
#ifdef DEBUG_SOLVEPROB
void solveProb::printResJac(int ioflag, int neq, const Array2D& Jac,
doublereal resid[], doublereal wtRes[],
@ -693,10 +630,7 @@ void solveProb::printResJac(int ioflag, int neq, const Array2D& Jac,
}
#endif
//================================================================================================
/*
* Optional printing at the start of the solveProb problem
*/
void solveProb::print_header(int ioflag, int ifunc, doublereal time_scale,
doublereal reltol,
doublereal netProdRate[])
@ -881,9 +815,8 @@ void solveProb::printIteration(int ioflag, doublereal damp, size_t label_d,
"------------------------------\n");
}
#endif
} /* printIteration */
}
//================================================================================================
void solveProb::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,
@ -971,7 +904,7 @@ void solveProb::printFinal(int ioflag, doublereal damp, size_t label_d, size_t l
}
#endif
}
//================================================================================================
#ifdef DEBUG_SOLVEPROB
void solveProb::
printIterationHeader(int ioflag, doublereal damp,doublereal inv_t, doublereal t_real,
@ -995,21 +928,19 @@ printIterationHeader(int ioflag, doublereal damp,doublereal inv_t, doublereal t_
}
}
#endif
//================================================================================================
void solveProb::setAtol(const doublereal atol[])
{
for (size_t k = 0; k < m_neq; k++, k++) {
m_atol[k] = atol[k];
}
}
//================================================================================================
void solveProb::setAtolConst(const doublereal atolconst)
{
for (size_t k = 0; k < m_neq; k++, k++) {
m_atol[k] = atolconst;
}
}
//================================================================================================
}