Doxygen update

This commit is contained in:
Harry Moffat 2011-01-11 17:34:31 +00:00
parent 3a63c938ea
commit f493f8f39a
5 changed files with 385 additions and 81 deletions

View file

@ -1,7 +1,8 @@
/**
* @file Integrator.h
*
* $Author$
*/
/* $Author$
* $Date$
* $Revision$
*

View file

@ -28,35 +28,37 @@
#include "mdp_allo.h"
#include <cfloat>
extern void print_line(const char *, int);
#include <vector>
#include <cstdio>
#include <cmath>
//@{
extern void print_line(const char *, int);
#ifndef MAX
#define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
#define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
#endif
//@}
using namespace std;
namespace Cantera {
//====================================================================================================================
//====================================================================================================================
//-----------------------------------------------------------
// Constants
//-----------------------------------------------------------
const doublereal DampFactor = 4;
//! 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
const int NDAMP = 7;
//====================================================================================================================
//-----------------------------------------------------------
// Static Functions
//-----------------------------------------------------------
//====================================================================================================================
//! Print a line of a single repeated character string
/*!
* @param str Character string
* @param n Iteration length
*/
static void print_line(const char *str, int n) {
for (int i = 0; i < n; i++) {
printf("%s", str);
@ -625,7 +627,8 @@ namespace Cantera {
* scaling has been implemented.
*/
int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr,
const doublereal * const ydot_curr, double* const delta_y, SquareMatrix& jac, int loglevel)
const doublereal * const ydot_curr, doublereal * const delta_y,
SquareMatrix& jac, int loglevel)
{
int irow;
@ -837,7 +840,6 @@ namespace Cantera {
return f_delta_bounds;
}
//====================================================================================================================
/*
*
* boundStep():
@ -865,7 +867,8 @@ namespace Cantera {
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*/
doublereal NonlinearSolver::boundStep(const doublereal * const y, const doublereal * const step0, const int loglevel) {
doublereal NonlinearSolver::boundStep(const doublereal * const y, const doublereal * const step0,
const int loglevel) {
int i, i_lower = -1;
doublereal fbound = 1.0, f_bounds = 1.0;
doublereal ff, y_new;
@ -1532,12 +1535,10 @@ namespace Cantera {
mdp::mdp_safe_free((void **) &imax);
}
//====================================================================================================================
/*
* subtractRD():
//! This routine subtracts two numbers for one another
/*!
* This routine subtracts 2 numbers. If the difference is less
* than 1.0E-14 times the magnitude of the smallest number,
* then diff returns an exact zero.
* than 1.0E-14 times the magnitude of the smallest number, then diff returns an exact zero.
* It also returns an exact zero if the difference is less than
* 1.0E-300.
*
@ -1546,8 +1547,12 @@ namespace Cantera {
* This routine is used in numerical differencing schemes in order
* to avoid roundoff errors resulting in creating Jacobian terms.
* Note: This is a slow routine. However, jacobian errors may cause
* loss of convergence. Therefore, in practice this routine
* has proved cost-effective.
* loss of convergence. Therefore, in practice this routine has proved cost-effective.
*
* @param a Value of a
* @param b value of b
*
* @return returns the difference between a and b
*/
static inline doublereal subtractRD(doublereal a, doublereal b) {
doublereal diff = a - b;

View file

@ -1,5 +1,5 @@
/**
* @file NonlinearSolve.h
* @file NonlinearSolver.h
* Class that calculates the solution to a nonlinear, dense, set
* of equations (see \ref numerics
* and class \link Cantera::NonlinearSolver NonlinearSolver\endlink).
@ -22,20 +22,55 @@
#include "ResidJacEval.h"
namespace Cantera {
// I think steady state is the only option I'm gunning for
//@{
/// @name Constant which determines the type of the nonlinear solve
/*!
* I think steady state is the only option I'm gunning for
*/
//! The nonlinear problem is part of a pseudo time dependent calculation (NOT TESTED)
#define NSOLN_TYPE_PSEUDO_TIME_DEPENDENT 2
//! The nonlinear problem is part of a time dependent calculation
#define NSOLN_TYPE_TIME_DEPENDENT 1
//! The nonlinear problem is part of a steady state calculation
#define NSOLN_TYPE_STEADY_STATE 0
//@}
//@{
/// @name Constant which determines the type of the Jacobian
//! The jacobian will be calculated from a numerical method
#define NSOLN_JAC_NUM 1
//! The jacobian is calculated from an analytical function
#define NSOLN_JAC_ANAL 2
//@}
//! Class that calculates the solution to a nonlinear system
/*!
* UNDER CONSTRUCTION - do not use!!!!!!!!!!!!!!!!!!!!!!!!!
* This is a small nonlinear solver that can solve highly nonlinear problems that
* must use a dense matrix to relax the system.
*
* Newton's method is used.
*
* Damping is used extensively when relaxing the system
*
*
*
* @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
*/
@ -77,7 +112,7 @@ namespace Cantera {
* The program always assumes that atol is specific
* to the solution component
*
* param y vector of the current solution values
* @param y vector of the current solution values
*/
void createSolnWeights(const doublereal * const y);
@ -155,9 +190,12 @@ namespace Cantera {
* recomputed. The row scales are recomputed here, after column
* scaling has been implemented.
*
* @param timeCurrent Current value of the time
* @param y_current Current value of the solution
* @param ydot_current Current value of the solution derivative.
* @param time_curr Current value of the time
* @param y_curr Current value of the solution
* @param ydot_curr Current value of the solution derivative.
* @param delta_y return value of the raw change in y
* @param jac Jacobian
* @param loglevel Log level
*
* @return Returns the result code from lapack. A zero means success. Anything
* else indicates a failure.
@ -178,7 +216,7 @@ namespace Cantera {
//! Set the delta Bounds magnitudes by hand
/*!
* @param deltaboundsMagnitudes
* @param deltaBoundsMagnitudes set the deltaBoundsMagnitude vector
*/
void setDeltaBoundsMagnitudes(const doublereal * const deltaBoundsMagnitudes);
@ -209,9 +247,13 @@ namespace Cantera {
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*
* @param y Current solution value of the old step
* @param step0 Proposed step change in the solution
* @param loglevel Log level
*
* @return Returns the damping factor determined by the bounds calculation
*/
doublereal boundStep(const double* const y, const double* const step0, const int loglevel);
doublereal boundStep(const doublereal * const y, const doublereal * const step0, const int loglevel);
//! Set bounds constraints for all variables in the problem
@ -242,10 +284,15 @@ namespace Cantera {
/*!
*
*
* @param J = Jacobian matrix to be filled in
* @param f = Right hand side. This routine returns the current
* @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.
* @param time_curr Current time
* @param CJ inverse of the value of deltaT
* @param y value of the solution vector
* @param ydot value of the time derivative of the solution vector
* @param num_newt_its Number of newton iterations
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
@ -259,6 +306,7 @@ namespace Cantera {
/*!
* @param timeCurrent Current value of the time
* @param ybase current value of the solution
* @param step0 Proposed step change in the solution
*
* @return Returns the norm of the value of the amount filtered
*/
@ -310,7 +358,16 @@ namespace Cantera {
* @param ydot0 Base value of the time derivative of teh
* solution
* @param step0 Initial step suggested.
* @param y1
* @param y1 Value of y1, the suggested solution after damping
* @param ydot1 Value of the time derivative of the solution at y1
* @param step1 Value of the step change from y0 to y1
* @param s1 norm of the step change in going from y0 to y1
* @param jac Jacobian
* @param loglevel Log level to be used
* @param writetitle Write a title line
* @param num_backtracks Number of backtracks taken
*
* @return returns an integer indicating what happened.
*/
int dampStep(const doublereal time_curr, const double* y0,
const doublereal *ydot0, const double* step0,
@ -331,6 +388,18 @@ namespace Cantera {
* equation system for now. Will make it more general later,
* if an application comes up.
*
* @param SolnType Solution type
* @param y_comm Initial value of the solution. On return this is the converged
* value of the solution
* @param ydot_comm Initial value of the solution derivative. On return this is the
* converged value of the solution derivative.
* @param CJ Inverse of the value of deltaT
* @param time_curr Current value of the time
* @param jac Matrix that will be used to store the jacobian
* @param num_newt_its Number of newton iterations taken
* @param num_linear_solves Number of linear solves taken
* @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
@ -389,6 +458,9 @@ namespace Cantera {
//! 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
*
* @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.
@ -429,6 +501,21 @@ namespace Cantera {
//! solution norms.
void calcSolnToResNormVector();
//! 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 setPrintLvl(int printLvl);
private:
@ -457,13 +544,17 @@ namespace Cantera {
//! Boolean indicating whether a manual delta steps have been input.
int m_manualDeltaStepSet;
//! Value of the delta step magnitudes
std::vector<doublereal> m_deltaStepMagnitudes;
//! Vector containing the current solution of the nonlinear solver
std::vector<doublereal> m_y_n;
//! Vector containing the solution at the previous time step
std::vector<doublereal> m_y_nm1;
//! New value of the solution time derivative
std::vector<doublereal> ydot_new;
//! Vector of column scaling factors
@ -483,9 +574,6 @@ namespace Cantera {
*/
std::vector<doublereal> m_rowWtScales;
//! Value of the residual for the nonlinear problem
std::vector<doublereal> m_resid;
@ -516,6 +604,7 @@ namespace Cantera {
//! Vector of the norm
doublereal m_normResidPoints[15];
//! Boolean indicating whether we should scale the residual
bool m_resid_scaled;
@ -523,7 +612,6 @@ namespace Cantera {
* INTERNAL BOUNDARY INFO FOR SOLUTIONS
*****************************************************************************************/
//! Bounds vector for each species
std::vector<doublereal> m_y_high_bounds;
@ -582,8 +670,10 @@ namespace Cantera {
*/
doublereal time_n;
//! Boolean indicating matrix conditioning
int m_matrixConditioning;
//! Order of the time step method = 1
int m_order;
//! value of the relative tolerance to use in solving the equation set
@ -592,8 +682,13 @@ namespace Cantera {
//! Base value of the absolute tolerance
doublereal atolBase_;
//! Vector containing the solution derivative at the previous time step
doublereal * m_ydot_nm1;
//! absolute tolerance in the solution unknown
/*!
* This is used to evaluating the weighting factor
*/
std::vector<doublereal> atolk_;
//! Determines the level of printing for each time step.
@ -620,10 +715,7 @@ namespace Cantera {
*/
static bool m_TurnOffTiming;
// Turn on or off printing of the Jacobian
/*!
*
*/
//! Turn on or off printing of the Jacobian
static bool s_print_NumJac;
};

View file

@ -134,50 +134,80 @@ namespace Cantera {
// Empty destructor
RootFind::~RootFind() {
}
//================================================================================================
double RootFind::delXNonzero(double x1) const {
double deltaX = 1.0E-14 * fabs(x1);
double delmin = DeltaXnorm_ * 1.0E-14;
//================================================================================================
// 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);
doublereal delmin = DeltaXnorm_ * 1.0E-14;
if (delmin > deltaX) {
return delmin;
}
return deltaX;
}
//================================================================================================
double RootFind::delXMeaningful(double x1) const {
double del = delXNonzero(x1);
// 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);
if (deltaXConverged_ > del) {
return deltaXConverged_;
}
return del;
}
//================================================================================================
double RootFind::deltaXControlled(double x2, double x1) const {
double sgnn = 1.0;
// Calcuated a controlled, nonzero delta between two numbers
/*
* The delta is designed to be greater than or equal to delXMeaningful(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 x2 first number
* @param x2 second number
*/
double RootFind::deltaXControlled(doublereal x2, doublereal x1) const {
doublereal sgnn = 1.0;
if (x1 > x2) {
sgnn = -1.0;
}
double deltaX = x2 - x1;
double x = fabs(x2) + fabs(x1);
double deltaXm = delXMeaningful(x);
doublereal deltaX = x2 - x1;
doublereal x = fabs(x2) + fabs(x1);
doublereal deltaXm = delXMeaningful(x);
if (fabs(deltaX) < deltaXm) {
deltaX = sgnn * deltaXm;
}
return deltaX;
}
//================================================================================================
bool RootFind::theSame(double x2, double x1) const {
double x = fabs(x2) + fabs(x1);
double deltaX = delXMeaningful(x);
//====================================================================================================================
// 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 delXMeaningful() of each other
*
* @param x2 First number
* @param x2 second number
*
* @return Returns a boolean indicating whether the two numbers are the same or not.
*/
bool RootFind::theSame(doublereal x2, doublereal x1) const {
doublereal x = fabs(x2) + fabs(x1);
doublereal deltaX = delXMeaningful(x);
if (fabs(x2 - x1) < deltaX) {
return true;
}
return false;
}
//================================================================================================
//====================================================================================================================
/*
* The following calculation is a line search method to find the root of a function
*
@ -313,8 +343,8 @@ namespace Cantera {
*/
foundStraddle = foundPosF && foundNegF;
if (foundStraddle) {
if (xPosF > xNegF) posStraddle = 1;
else posStraddle = 0 ;
if (xPosF > xNegF) posStraddle = 1;
else posStraddle = 0 ;
}
bool doQuad = false;
bool useNextStrat = false;
@ -332,7 +362,7 @@ namespace Cantera {
printf(" RootFind: we are here x2 = %g x1 = %g\n", x2, x1);
}
#endif
double delXtmp = deltaXControlled(x2, x1);
doublereal delXtmp = deltaXControlled(x2, x1);
slope = (f2 - f1) / delXtmp;
if (fabs(slope) <= 1.0E-100) {
if (printLvl >= 2) {
@ -715,7 +745,7 @@ namespace Cantera {
return retn;
}
//================================================================================================
//====================================================================================================================
doublereal RootFind::func(doublereal x) {
doublereal r;
#ifdef DEBUG_MODE
@ -727,18 +757,52 @@ namespace Cantera {
#endif
return (r - m_funcTargetValue);
}
//================================================================================================
//====================================================================================================================
// 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 rtol, doublereal atol)
{
m_atol = atol;
m_rtol = rtol;
}
//================================================================================================
//====================================================================================================================
// 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) {
@ -746,7 +810,18 @@ namespace Cantera {
}
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) {
@ -754,10 +829,16 @@ namespace Cantera {
}
FuncIsGenerallyDecreasing_ = value;
}
//================================================================================================
//====================================================================================================================
// Set the minimum value of deltaX
/*
* This sets the value of deltaXNorm_
*
* @param deltaXNorm
*/
void RootFind::setDeltaX(doublereal deltaXNorm)
{
DeltaXnorm_ = deltaXNorm;
}
//================================================================================================
//====================================================================================================================
}

View file

@ -24,9 +24,16 @@
namespace Cantera {
//@{
/// @name Constant which determines the return integer from the routine
//! This means that the root solver was a success
#define ROOTFIND_SUCCESS 0
//! This means that the root solver failed to achieve convergence
#define ROOTFIND_FAILEDCONVERGENCE -1
//! This means that the input to the root solver was defective
#define ROOTFIND_BADINPUT -2
//@{
//! Root finder for 1D problems
/*!
@ -40,6 +47,8 @@ namespace Cantera {
//! Constructor for the object
/*!
*
* @param resid Pointer to the residual function to be used to calculate f(x)
*/
RootFind(ResidEval* resid);
@ -49,19 +58,57 @@ namespace Cantera {
private:
//! Unimplemented private copy constructor
/*!
* @param right object to be copied
*/
RootFind(const RootFind &right);
//! Unimplemented private assignment operator
/*!
* @param right object to be copied
*/
RootFind& operator=(const RootFind &right);
//! 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 delXNonzero(doublereal x1) const;
double delXNonzero(double x1) const;
double delXMeaningful(double x1) const;
//! 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 delXMeaningful(doublereal x1) const;
double deltaXControlled(double x2, double x1) const;
//! Calcuated a controlled, nonzero delta between two numbers
/*!
* The delta is designed to be greater than or equal to delXMeaningful(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 x2 first number
* @param x1 second number
*/
doublereal deltaXControlled(doublereal x2, doublereal x1) const;
bool theSame(double x2, double x1) const;
//! 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 delXMeaningful() of each other
*
* @param x1 First number
* @param x2 second number
*
* @return Returns a boolean indicating whether the two numbers are the same or not.
*/
bool theSame(doublereal x2, doublereal x1) const;
public:
@ -105,24 +152,102 @@ namespace Cantera {
*/
doublereal func(doublereal x);
//! 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 setTol(doublereal rtol, doublereal atol);
//! 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 setPrintLvl(int 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 setFuncIsGenerallyIncreasing(bool 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 setFuncIsGenerallyDecreasing(bool value);
//! Set the minimum value of deltaX
/*!
* This sets the value of deltaXNorm_
*
* @param deltaXNorm
*/
void setDeltaX(doublereal deltaXNorm);
public:
//! Pointer to the residual function evaluator
ResidEval *m_residFunc;
//! Target value for the function. We seek the value of f that is equal to this value
doublereal m_funcTargetValue;
//! Absolute tolerance for the value of f
doublereal m_atol;
//! Relative tolerance for the value of f
doublereal m_rtol;
//! Maximum number of step sizes
doublereal m_maxstep;
protected:
//! Print level
int printLvl;
//! Delta X norm. This is the minimum value of deltaX that will be used by the program
doublereal DeltaXnorm_;
//! Boolean indicating whether the function is an increasing with x
bool FuncIsGenerallyIncreasing_;
//! Boolean indicating whether the function is decreasing with x
bool FuncIsGenerallyDecreasing_;
//! Value of delta X that is needed for convergence
/*!
* X will be considered as converged if we are within deltaXConverged_ of the solution
* The default is zero.
*/
doublereal deltaXConverged_;
};