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

Added return variables to a lot of Residual evaluation functions.

Browse source 
These are now used by the nonlinear solver to bail if the residual
says to bail.
This commit is contained in:
Harry Moffat 2011-01-04 22:04:59 +00:00
parent 5395e214b8
commit b1750d579a
7 changed files with 258 additions and 81 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

50
Cantera/src/numerics/NonlinearSolver.cpp
View file

@ -487,13 +487,18 @@ namespace Cantera {
* @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
*/
void NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr,
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)
{
m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, DATA_PTR(m_resid), evalType);
int retn = m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, DATA_PTR(m_resid), evalType);
m_nfe++;
m_resid_scaled = false;
return retn;
}
//====================================================================================================================
// Compute the undamped Newton step
@ -994,9 +999,15 @@ namespace Cantera {
* -> m_resid[] contains the result of the residual calculation
*/
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
doResidualCalc(time_curr, solnType_, y1, ydot1, Base_LaggedSolutionComponents);
info = doResidualCalc(time_curr, solnType_, y1, ydot1, Base_LaggedSolutionComponents);
} else {
doResidualCalc(time_curr, solnType_, y1, ydot0, Base_LaggedSolutionComponents);
info = doResidualCalc(time_curr, solnType_, y1, ydot0, Base_LaggedSolutionComponents);
}
if (info != 1) {
if (loglevel > 0) {
printf("\t\t\tdampStep: current trial step and damping led to Residual Calc ERROR %d. Bailing\n", info);
}
return -1;
}
m_normResidTrial = residErrorNorm(DATA_PTR(m_resid));
@ -1143,8 +1154,10 @@ namespace Cantera {
*
* 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.
*
* if an application comes up.
*
* @return A positive value indicates a successful convergence
* -1 Failed convergence
*/
int NonlinearSolver::solve_nonlinear_problem(int SolnType, double* y_comm,
double* ydot_comm, doublereal CJ,
@ -1261,8 +1274,14 @@ namespace Cantera {
doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(m_y_n), DATA_PTR(ydot_curr));
info = doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(m_y_n), DATA_PTR(ydot_curr));
if (info != 1) {
if (m_print_flag > 0) {
printf("\t\t\tsolve_nonlinear_problem(): Residual Calc ERROR %d. Bailing\n", info);
}
m = -1;
goto done;
}
/*
* Scale the matrix and the rhs, if they aren't already scaled
@ -1340,8 +1359,14 @@ namespace Cantera {
}
doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(y_new), DATA_PTR(ydot_new));
info = doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(y_new), DATA_PTR(ydot_new));
if (info != 1) {
if (m_print_flag > 0) {
printf("\t\t\tdampStep: current trial step and damping led to Residual Calc ERROR %d. Bailing\n", info);
}
m = -1;
goto done;
}
if (m_print_flag > 3) {
residErrorNorm(DATA_PTR(m_resid), "Resulting Residual Norm", 10, DATA_PTR(y_new));
@ -1884,5 +1909,10 @@ namespace Cantera {
}
//=====================================================================================================================
void NonlinearSolver::setPrintLvl( int printLvl)
{
m_print_flag = printLvl;
}
//=====================================================================================================================
}

52
Cantera/src/numerics/NonlinearSolver.h
View file

@ -96,6 +96,8 @@ namespace Cantera {
* @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.
*
* @return Returns the L2 norm of the delta
*/
doublereal solnErrorNorm(const doublereal * const delta_y, const char * title = 0, int printLargest = 0,
const doublereal dampFactor = 1.0);
@ -113,26 +115,13 @@ namespace Cantera {
* @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,
const doublereal * const y = 0);
//! Compute the current Residual
/*!
* Compute the time dependent residual of
* the set of equations.
*/
// void doTDResidualCalc(const double time_curr, const int typeCalc,
// const double * const y_curr, const double * const ydot_curr, int loglevel);
//! Compute the current Residual
/*!
* Compute the steady state residual of
* the set of equations.
*/
// void doSteadyResidualCalc(const double time_curr, const int typeCalc,
// const double * const y_curr, int loglevel);
//! Compute the current residual
/*!
* The current value of the residual is storred in the internal work array m_resid.
@ -143,10 +132,14 @@ namespace Cantera {
* @param ydot_curr Current value of the time derivative of the solution vector
* @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
*/
void doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr,
const doublereal * const ydot_curr,
const ResidEval_Type_Enum evalType = Base_ResidEval);
int doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr,
const doublereal * const ydot_curr,
const ResidEval_Type_Enum evalType = Base_ResidEval);
//! Compute the undamped Newton step
/*!
@ -166,7 +159,7 @@ namespace Cantera {
* @param y_current Current value of the solution
* @param ydot_current Current value of the solution derivative.
*
* @return returns the result code from lapack. A zero means success. Anything
* @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,
@ -215,9 +208,10 @@ namespace Cantera {
* factor of 2
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*
* @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 double* const y, const double* const step0, const int loglevel);
//! Set bounds constraints for all variables in the problem
@ -229,13 +223,13 @@ namespace Cantera {
void setBoundsConstraints(const doublereal * const y_low_bounds,
const doublereal * const y_high_bounds);
//! return an editable vector of the low bounds constraints
//! Return an editable vector of the low bounds constraints
std::vector<double> & lowBoundsConstraintVector();
//! return an editable vector of the high bounds constraints
//! Return an editable vector of the high bounds constraints
std::vector<double> & highBoundsConstraintVector();
//! Internal function to calculate the time derivative at the new step
//! Internal function to calculate the time derivative at the new step
/*!
* @param order of the BDF method
* @param y_curr current value of the solution
@ -256,13 +250,11 @@ namespace Cantera {
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* 0 Means an unsuccessful operation
*
*/
int beuler_jac(SquareMatrix &J, doublereal * const f,
doublereal time_curr, doublereal CJ, doublereal * const y,
doublereal * const ydot, int num_newt_its);
//! Apply a filtering process to the step
/*!
* @param timeCurrent Current value of the time
@ -290,9 +282,9 @@ namespace Cantera {
* 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
* factor of 1.5
* Maximum decrease in variable in any one newton iteration:
* factor of 2
* factor of 2
*
* @param y Initial value of the solution vector
* @param step0 initial proposed step size
@ -437,6 +429,8 @@ namespace Cantera {
//! solution norms.
void calcSolnToResNormVector();
void setPrintLvl(int printLvl);
private:
//! Pointer to the residual and jacobian evaluator for the

5
Cantera/src/numerics/ResidEval.h
View file

@ -101,10 +101,13 @@ namespace Cantera {
/**
* 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
*/
virtual void getInitialConditions(const doublereal t0, doublereal * const y,
virtual int getInitialConditions(const doublereal t0, doublereal * const y,
doublereal * const ydot) {
throw CanteraError("ResidEval::GetInitialConditions()", "base class called");
return 1;
}
//! Return the number of equations in the equation system

15
Cantera/src/numerics/ResidJacEval.cpp
View file

@ -96,7 +96,7 @@ namespace Cantera {
* @param y Solution vector (output)
* @param ydot Rate of change of solution vector. (output)
*/
void ResidJacEval::
int ResidJacEval::
getInitialConditions(doublereal t0, doublereal * const y, doublereal * const ydot) {
for (int i = 0; i < neq_; i++) {
y[i] = 0.0;
@ -106,6 +106,7 @@ namespace Cantera {
ydot[i] = 0.0;
}
}
return 1;
}
//====================================================================================================================
// This function may be used to create output at various points in the execution of an application.
@ -154,10 +155,11 @@ namespace Cantera {
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
*/
void ResidJacEval::
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
@ -279,9 +281,10 @@ namespace Cantera {
* @param nrows offsets for the matrix
* @param rhs residual vector. This also needs to be lhs multiplied by M
*/
void ResidJacEval::
int ResidJacEval::
matrixConditioning(doublereal * const matrix, const int nrows, doublereal * const rhs)
{
return 1;
}
//====================================================================================================================
// Evaluate the residual function
@ -297,12 +300,13 @@ namespace Cantera {
* differenced or that the residual doesn't take this issue into account)
* @param delta_x Value of the delta used in the numerical differencing
*/
void ResidJacEval::
int ResidJacEval::
evalResidNJ(const doublereal t, const doublereal deltaT, const doublereal * y,
const doublereal * ydot, doublereal * const resid, const ResidEval_Type_Enum evalType,
const int id_x, const doublereal delta_x)
{
throw CanteraError("ResidJacEval::evalResidNJ()", "Not implemented\n");
return 1;
}
//====================================================================================================================
// Calculate an analytical jacobian and the residual at the current time and values.
@ -316,7 +320,7 @@ namespace Cantera {
* @param J Reference to the SquareMatrix object to be calculated (output)
* @param resid Value of the residual that is computed (output)
*/
void ResidJacEval::
int ResidJacEval::
evalJacobian(const doublereal t, const doublereal delta_t,
const doublereal * const y,
const doublereal * const ydot,
@ -324,6 +328,7 @@ namespace Cantera {
doublereal * const resid)
{
throw CanteraError("ResidJacEval::evalJacobian()", "Not implemented\n");
return 1;
}
//====================================================================================================================

38
Cantera/src/numerics/ResidJacEval.h
View file

@ -112,8 +112,12 @@ namespace Cantera {
* 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
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
virtual void evalResidNJ(const doublereal t, const doublereal delta_t,
virtual int evalResidNJ(const doublereal t, const doublereal delta_t,
const doublereal * const y,
const doublereal * const ydot,
doublereal * const resid,
@ -129,8 +133,12 @@ namespace Cantera {
* @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 void getInitialConditions(const doublereal t0, doublereal * const y, doublereal * const ydot);
virtual int getInitialConditions(const doublereal t0, doublereal * const y, doublereal * const ydot);
//! Filter the solution predictions
/*!
@ -174,8 +182,12 @@ namespace Cantera {
* @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)
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
virtual void evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal * const y,
virtual int evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal * const y,
const doublereal * const ydot);
//! Evalulate any stopping criteria other than a final time limit
@ -281,8 +293,12 @@ namespace Cantera {
* @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
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
virtual void matrixConditioning(doublereal * const matrix, const int nrows,
virtual int matrixConditioning(doublereal * const matrix, const int nrows,
doublereal * const rhs);
//! Calculate an analytical jacobian and the residual at the current time and values.
@ -295,12 +311,16 @@ namespace Cantera {
* @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)
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
virtual void evalJacobian(const doublereal t, const doublereal delta_t,
const doublereal* const y,
const doublereal* const ydot,
SquareMatrix &J,
doublereal * const resid);
virtual int evalJacobian(const doublereal t, const doublereal delta_t,
const doublereal* const y,
const doublereal* const ydot,
SquareMatrix &J,
doublereal * const resid);
protected:

125
Cantera/src/numerics/RootFind.cpp
View file

@ -125,7 +125,8 @@ namespace Cantera {
printLvl(0),
DeltaXnorm_(0.01),
FuncIsGenerallyIncreasing_(false),
FuncIsGenerallyDecreasing_(false)
FuncIsGenerallyDecreasing_(false),
deltaXConverged_(0.0)
{
}
@ -133,6 +134,49 @@ 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;
if (delmin > deltaX) {
return delmin;
}
return deltaX;
}
//================================================================================================
double RootFind::delXMeaningful(double x1) const {
double del = delXNonzero(x1);
if (deltaXConverged_ > del) {
return deltaXConverged_;
}
return del;
}
//================================================================================================
double RootFind::deltaXControlled(double x2, double x1) const {
double sgnn = 1.0;
if (x1 > x2) {
sgnn = -1.0;
}
double deltaX = x2 - x1;
double x = fabs(x2) + fabs(x1);
double 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);
if (fabs(x2 - x1) < deltaX) {
return true;
}
return false;
}
//================================================================================================
/*
* The following calculation is a line search method to find the root of a function
@ -162,6 +206,7 @@ namespace Cantera {
FILE *fp = 0;
#endif
doublereal x1, x2, xnew, f1, f2, fnew, slope;
doublereal deltaX1 = 0.0, deltaX2 = 0.0, deltaXnew = 0.0;
int its = 0;
int posStraddle = 0;
int retn = 0;
@ -173,6 +218,7 @@ namespace Cantera {
doublereal fnorm; /* A valid norm for the making the function value dimensionless */
doublereal c[9], f[3], xn1, xn2, x0 = 0.0, f0 = 0.0, root, theta, xquad, xDelMin;
doublereal CR0, CR1, CR2, CRnew, CRdenom;
doublereal sgn;
callNum++;
#ifdef DEBUG_MODE
@ -230,6 +276,7 @@ namespace Cantera {
x2 = x1 - (xmax - xmin) / 100.;
}
deltaX2 = x2 - x1;
f2 = func(x2);
#ifdef DEBUG_MODE
if (printLvl >= 3) {
@ -277,7 +324,13 @@ namespace Cantera {
* Find an estimate of the next point, xnew, to try based on
* a linear approximation from the last two points.
*/
slope = (f2 - f1) / (x2 - x1);
#ifdef DEBUG_HKM
if (fabs(x2 - x1) < 1.0E-14) {
printf(" RootFind: we are here x2 = %g x1 = %g\n", x2, x1);
}
#endif
double delXtmp = deltaXControlled(x2, x1);
slope = (f2 - f1) / delXtmp;
if (fabs(slope) <= 1.0E-100) {
if (printLvl >= 2) {
writelogf("%s functions evals produced the same result, %g, at %g and %g\n",
@ -299,6 +352,7 @@ namespace Cantera {
fprintf(fp, " | xlin = %-11.5E", xnew);
}
#endif
deltaXnew = xnew - x2;
/*
* If the suggested step size is too big, throw out step
*/
@ -306,6 +360,10 @@ namespace Cantera {
if (fabs(xnew - x2) > 3.0 * DeltaXnorm_) {
useNextStrat = true;
}
if (fabs(deltaXnew) < fabs(deltaX2)) {
deltaXnew = DSIGN(deltaXnew) * 1.1 * fabs(deltaX2);
xnew = deltaXnew + x2;
}
}
if (useNextStrat) {
if (f2 < 0.0) {
@ -442,26 +500,7 @@ namespace Cantera {
#endif
}
}
/*
* Guard against going above xmax or below xmin
*/
if (xnew > xmax) {
xnew = x2 + (xmax - x2) / 2.0;
#ifdef DEBUG_MODE
if (printLvl >= 3) {
fprintf(fp, " | xlimitmax = %-11.5E", xnew);
}
#endif
}
if (xnew < xmin) {
xnew = x2 + (x2 - xmin) / 2.0;
#ifdef DEBUG_MODE
if (printLvl >= 3) {
fprintf(fp, " | xlimitmin = %-11.5E", xnew);
}
#endif
}
if (foundStraddle) {
#ifdef DEBUG_MODE
slope = xnew;
@ -507,7 +546,42 @@ namespace Cantera {
}
#endif
}
/*
* Enforce a minimum stepsize if we haven't found a straddle.
*/
deltaXnew = xnew - x2;
if (fabs(deltaXnew) < 1.2 * delXMeaningful(xnew)) {
if (!foundStraddle) {
sgn = 1.0;
if (x2 < xnew) {
sgn = -1.0;
}
deltaXnew = 1.2 * delXMeaningful(xnew) * sgn;
xnew = x2 + deltaXnew;
}
}
/*
* Guard against going above xmax or below xmin
*/
if (xnew > xmax) {
xnew = x2 + (xmax - x2) / 2.0;
#ifdef DEBUG_MODE
if (printLvl >= 3) {
fprintf(fp, " | xlimitmax = %-11.5E", xnew);
}
#endif
}
if (xnew < xmin) {
xnew = x2 + (x2 - xmin) / 2.0;
#ifdef DEBUG_MODE
if (printLvl >= 3) {
fprintf(fp, " | xlimitmin = %-11.5E", xnew);
}
#endif
}
fnew = func(xnew);
CRdenom = MAX(fabs(fnew), MAX(fabs(f2), MAX(fabs(f1), fnorm)));
CRnew = sqrt(fabs(fnew) / CRdenom);
@ -563,6 +637,8 @@ namespace Cantera {
CR1 = CR2;
x2 = xnew;
f2 = fnew;
deltaX1 = deltaX2;
deltaX2 = deltaXnew;
CR2 = CRnew;
if (fabs(fnew / fnorm) < m_rtol) {
converged = 1;
@ -576,10 +652,9 @@ namespace Cantera {
retn = ROOTFIND_FAILEDCONVERGENCE;
converged = true;
}
if (fabs(x2 - x1) / denom < 1.0E-13) {
if (theSame(x2, x1)) {
converged = true;
}
}
its++;
} while (! converged && its < itmax);

54
Cantera/src/numerics/RootFind.h
View file

@ -28,6 +28,12 @@ namespace Cantera {
#define ROOTFIND_FAILEDCONVERGENCE -1
#define ROOTFIND_BADINPUT -2
//! Root finder for 1D problems
/*!
*
*
*
*/
class RootFind {
public:
@ -48,10 +54,53 @@ namespace Cantera {
//! Unimplemented private assignment operator
RootFind& operator=(const RootFind &right);
public:
int solve(doublereal xmin, doublereal xmax, int itmax, doublereal funcTargetValue, doublereal *xbest) ;
double delXNonzero(double x1) const;
double delXMeaningful(double x1) const;
double deltaXControlled(double x2, double x1) const;
bool theSame(double x2, double x1) const;
public:
//! Using a line search method, find the root of a 1D function
/*!
* This routine solves the following equation.
*
* \f[
* R(x) = f(x) - f_o = 0
* \f]
*
* @param xmin Minimum value of x to be used.
* @param xmax Maximum value of x to be used
* @param itmax maximum number of iterations. Usually, it can be less than 50.
* @param funcTargetValue Value of \f$ f_o \f$ in the equation.
* @param xbest Returns the x that satisfies the function
* On input, xbest should contain the best estimate of the solution.
* An attempt to find the solution near xbest is made.
*
* @return:
* 0 = ROOTFIND_SUCCESS Found function
* -1 = ROOTFIND_FAILEDCONVERGENCE Failed to find the answer
* -2 = ROOTFIND_BADINPUT Bad input was detected
*/
int solve(doublereal xmin, doublereal xmax, int itmax, doublereal funcTargetValue, doublereal *xbest);
//! Return the function value
/*!
* This routine evaluates the following equation.
*
* \f[
* R(x) = f(x) - f_o = 0
* \f]
*
* @param x Value of the independent variable
*
* @return The routine returns the value of \f$ R(x) \f$
*/
doublereal func(doublereal x);
void setTol(doublereal rtol, doublereal atol);
@ -72,6 +121,7 @@ namespace Cantera {
doublereal DeltaXnorm_;
bool FuncIsGenerallyIncreasing_;
bool FuncIsGenerallyDecreasing_;
doublereal deltaXConverged_;
};
}