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

Added hooks for rtol and max newton iterations

Browse source
This commit is contained in:
Harry Moffat 2010-10-14 22:34:59 +00:00
parent 97f326a628
commit fad1f9f7ff
2 changed files with 157 additions and 124 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -50,7 +50,7 @@ namespace Cantera {
// Constants
//-----------------------------------------------------------
const double DampFactor = 4;
const doublereal DampFactor = 4;
const int NDAMP = 7;
//====================================================================================================================
//-----------------------------------------------------------
@ -103,6 +103,7 @@ namespace Cantera {
m_numTotalNewtIts(0),
m_min_newt_its(0),
filterNewstep(0),
maxNewtIts_(50),
m_jacFormMethod(NSOLN_JAC_NUM),
m_nJacEval(0),
time_n(0.0),
@ -171,6 +172,7 @@ namespace Cantera {
m_numTotalNewtIts(0),
m_min_newt_its(0),
filterNewstep(0),
maxNewtIts_(50),
m_jacFormMethod(NSOLN_JAC_NUM),
m_nJacEval(0),
time_n(0.0),
@ -228,6 +230,7 @@ namespace Cantera {
m_numTotalNewtIts = right.m_numTotalNewtIts;
m_min_newt_its = right.m_min_newt_its;
filterNewstep = right.filterNewstep;
maxNewtIts_ = right.maxNewtIts_;
m_jacFormMethod = right.m_jacFormMethod;
m_nJacEval = right.m_nJacEval;
time_n = right.time_n;
@ -252,7 +255,7 @@ namespace Cantera {
*
* param y vector of the current solution values
*/
void NonlinearSolver::createSolnWeights(const double * const y) {
void NonlinearSolver::createSolnWeights(const doublereal * const y) {
for (int i = 0; i < neq_; i++) {
m_ewt[i] = rtol_ * fabs(y[i]) + atolk_[i];
}
@ -264,8 +267,8 @@ namespace Cantera {
* @param y_low_bounds Vector of lower bounds
* @param y_high_bounds Vector of high bounds
*/
void NonlinearSolver::setBoundsConstraints(const double * const y_low_bounds,
const double * const y_high_bounds) {
void NonlinearSolver::setBoundsConstraints(const doublereal * const y_low_bounds,
const doublereal * const y_high_bounds) {
for (int i = 0; i < neq_; i++) {
m_y_low_bounds[i] = y_low_bounds[i];
m_y_high_bounds[i] = y_high_bounds[i];
@ -287,11 +290,11 @@ namespace Cantera {
* @param dampFactor Current value of the damping factor. Defaults to 1.
* only used for printout out a table.
*/
double NonlinearSolver::solnErrorNorm(const double * const delta_y, const char * title, int printLargest,
const double dampFactor)
doublereal NonlinearSolver::solnErrorNorm(const doublereal * const delta_y, const char * title, int printLargest,
const doublereal dampFactor)
{
int i;
double sum_norm = 0.0, error;
doublereal sum_norm = 0.0, error;
for (i = 0; i < neq_; i++) {
error = delta_y[i] / m_ewt[i];
sum_norm += (error * error);
@ -322,7 +325,7 @@ namespace Cantera {
}
printf(" = %-11.4E\n", sum_norm);
double dmax1, normContrib;
doublereal dmax1, normContrib;
int j;
int *imax = mdp::mdp_alloc_int_1(num_entries, -1);
printf("\t\t Printout of Largest Contributors:\n");
@ -372,11 +375,11 @@ namespace Cantera {
* if true, then a table of the largest values is printed
* out to standard output.
*/
double NonlinearSolver::residErrorNorm(const double * const resid, const char * title, const int printLargest,
const double * const y)
doublereal NonlinearSolver::residErrorNorm(const doublereal * const resid, const char * title, const int printLargest,
const doublereal * const y)
{
int i;
double sum_norm = 0.0, error;
doublereal sum_norm = 0.0, error;
for (i = 0; i < neq_; i++) {
error = resid[i] / m_residWts[i];
sum_norm += (error * error);
@ -384,7 +387,7 @@ namespace Cantera {
sum_norm = sqrt(sum_norm / neq_);
if (printLargest) {
const int num_entries = printLargest;
double dmax1, normContrib;
doublereal dmax1, normContrib;
int j;
int *imax = mdp::mdp_alloc_int_1(num_entries, -1);
@ -461,8 +464,8 @@ namespace Cantera {
* @param y_curr Current value of the solution vector
* @param ydot_curr Current value of the time derivative of the solution vector
*/
void NonlinearSolver::doResidualCalc(const double time_curr, const int typeCalc, const double * const y_curr,
const double * const ydot_curr)
void NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr,
const doublereal * const ydot_curr)
{
m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, DATA_PTR(m_resid), Base_ResidEval);
m_nfe++;
@ -482,7 +485,7 @@ namespace Cantera {
* recomputed. The row scales are recomputed here, after column
* scaling has been implemented.
*/
void NonlinearSolver::scaleMatrix(SquareMatrix& jac, double* y_comm, double* ydot_comm, double time_curr)
void NonlinearSolver::scaleMatrix(SquareMatrix& jac, double* y_comm, double* ydot_comm, doublereal time_curr)
{
int irow, jcol;
@ -503,7 +506,7 @@ namespace Cantera {
/*
* Scale the new Jacobian
*/
double *jptr = &(*(jac.begin()));
doublereal *jptr = &(*(jac.begin()));
for (jcol = 0; jcol < neq_; jcol++) {
for (irow = 0; irow < neq_; irow++) {
*jptr *= m_colScales[jcol];
@ -523,9 +526,9 @@ namespace Cantera {
/*
* Ok, this is ugly. jac.begin() returns an vector<double> iterator
* to the first data location.
* Then &(*()) reverts it to a double *.
* Then &(*()) reverts it to a doublereal *.
*/
double *jptr = &(*(jac.begin()));
doublereal *jptr = &(*(jac.begin()));
for (irow = 0; irow < neq_; irow++) {
m_rowScales[irow] = 0.0;
m_rowWtScales[irow] = 0.0;
@ -578,8 +581,8 @@ namespace Cantera {
* recomputed. The row scales are recomputed here, after column
* scaling has been implemented.
*/
void NonlinearSolver::doNewtonSolve(const double time_curr, const double * const y_curr,
const double * const ydot_curr, double* const delta_y, SquareMatrix& jac, int loglevel)
void NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr,
const doublereal * const ydot_curr, double* const delta_y, SquareMatrix& jac, int loglevel)
{
int irow;
@ -618,9 +621,9 @@ namespace Cantera {
if (printJacContributions) {
for (int iNum = 0; iNum < numRows; iNum++) {
if (iNum > 0) focusRow++;
double dsum = 0.0;
doublereal dsum = 0.0;
vector_fp& Jdata = jacBack.data();
double dRow = Jdata[neq_ * focusRow + focusRow];
doublereal dRow = Jdata[neq_ * focusRow + focusRow];
printf("\n Details on delta_Y for row %d \n", focusRow);
printf(" Value before = %15.5e, delta = %15.5e,"
"value after = %15.5e\n", y_curr[focusRow],
@ -639,8 +642,8 @@ namespace Cantera {
dsum += RRow[iNum] / dRow;
for (int ii = 0; ii < neq_; ii++) {
if (ii != focusRow) {
double aij = Jdata[neq_ * ii + focusRow];
double contrib = aij * delta_y[ii] * (-1.0) / dRow;
doublereal aij = Jdata[neq_ * ii + focusRow];
doublereal contrib = aij * delta_y[ii] * (-1.0) / dRow;
dsum += contrib;
if (fabs(contrib) > Pcutoff) {
printf("%6d %15.5e %15.5e %15.5e\n", ii,
@ -668,7 +671,7 @@ namespace Cantera {
}
}
//====================================================================================================================
void NonlinearSolver::setDeltaBoundsMagnitudes(const double * const deltaBoundsMagnitudes)
void NonlinearSolver::setDeltaBoundsMagnitudes(const doublereal * const deltaBoundsMagnitudes)
{
for (int i = 0; i < neq_; i++) {
@ -696,18 +699,18 @@ namespace Cantera {
* @return returns the damping factor
*/
double
NonlinearSolver::deltaBoundStep(const double * const y, const double * const step0, const int loglevel) {
NonlinearSolver::deltaBoundStep(const doublereal * const y, const doublereal * const step0, const int loglevel) {
int i_fbounds = 0;
int ifbd = 0;
int i_fbd = 0;
double sameSign = 0.0;
double ff;
double f_delta_bounds = 1.0;
double ff_alt;
doublereal sameSign = 0.0;
doublereal ff;
doublereal f_delta_bounds = 1.0;
doublereal ff_alt;
for (int i = 0; i < neq_; i++) {
double y_new = y[i] + step0[i];
doublereal y_new = y[i] + step0[i];
sameSign = y_new * y[i];
/*
@ -818,10 +821,10 @@ namespace Cantera {
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*/
double NonlinearSolver::boundStep(const double * const y, const double * const step0, const int loglevel) {
doublereal NonlinearSolver::boundStep(const doublereal * const y, const doublereal * const step0, const int loglevel) {
int i, i_lower = -1;
double fbound = 1.0, f_bounds = 1.0;
double ff, y_new;
doublereal fbound = 1.0, f_bounds = 1.0;
doublereal ff, y_new;
for (i = 0; i < neq_; i++) {
y_new = y[i] + step0[i];
@ -830,7 +833,7 @@ namespace Cantera {
*/
if (step0[i] < 0.0) {
if (y_new < (y[i] + 0.8 * (m_y_low_bounds[i] - y[i]))) {
double legalDelta = 0.8*(m_y_low_bounds[i] - y[i]);
doublereal legalDelta = 0.8*(m_y_low_bounds[i] - y[i]);
ff = legalDelta / step0[i];
if (ff < f_bounds) {
f_bounds = ff;
@ -843,7 +846,7 @@ namespace Cantera {
*/
if (step0[i] > 0.0) {
if (y_new > (y[i] + 0.8 * (m_y_high_bounds[i] - y[i]))) {
double legalDelta = 0.8*(m_y_high_bounds[i] - y[i]);
doublereal legalDelta = 0.8*(m_y_high_bounds[i] - y[i]);
ff = legalDelta / step0[i];
if (ff < f_bounds) {
f_bounds = ff;
@ -864,7 +867,7 @@ namespace Cantera {
}
}
double f_delta_bounds = deltaBoundStep(y, step0, loglevel);
doublereal f_delta_bounds = deltaBoundStep(y, step0, loglevel);
fbound = MIN(f_bounds, f_delta_bounds);
return fbound;
@ -886,23 +889,23 @@ namespace Cantera {
* 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 deltaSoln is not produced however. s1 is estimated.
* 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 deltaSoln is not produced however. s1 is estimated.
* A predicted deltaSoln1 is not produced however. s1 is estimated.
* 0 Uncertain Success: s1 is about the same as s0
* -2 Unsuccessful step.
*/
int NonlinearSolver::dampStep(const double time_curr, const double* y0,
const double *ydot0, const double* step0,
int NonlinearSolver::dampStep(const doublereal time_curr, const double* y0,
const doublereal *ydot0, const double* step0,
double* const y1, double* const ydot1, double* step1,
double& s1, SquareMatrix& jac,
int& loglevel, bool writetitle,
int& num_backtracks) {
double& s1, SquareMatrix& jac, int& loglevel, bool writetitle,
int& num_backtracks)
{
int retnTrial = -2;
// Compute the weighted norm of the undamped step size step0
double s0 = solnErrorNorm(step0);
doublereal s0 = solnErrorNorm(step0);
// Compute the multiplier to keep all components in bounds
// A value of one indicates that there is no limitation
@ -927,7 +930,7 @@ namespace Cantera {
// damping coefficient starts at 1.0
m_dampRes = 1.0;
int j, m;
double ff = m_dampBound;
doublereal ff = m_dampBound;
num_backtracks = 0;
for (m = 0; m < NDAMP; m++) {
@ -947,6 +950,10 @@ namespace Cantera {
} else {
}
/*
* Calculate the residual that would result if y1[] were the new solution vector
* -> m_resid[] contains the result of the residual calculation
*/
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
doResidualCalc(time_curr, solnType_, y1, ydot1);
} else {
@ -988,8 +995,8 @@ namespace Cantera {
break;
}
// compute the next undamped step, step1[], that would result
// if y1[] were accepted.
// Compute the next undamped step, step1[], that would result if y1[] were accepted.
// We now have two steps that we have calculated step0[] and step1[]
if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
doNewtonSolve(time_curr, y1, ydot1, step1, jac, loglevel);
} else {
@ -1001,9 +1008,9 @@ namespace Cantera {
// write log information
if (loglevel > 3) {
print_solnDelta_norm_contrib((const double *) step0,
print_solnDelta_norm_contrib((const doublereal *) step0,
"DeltaSoln",
(const double *) step1,
(const doublereal *) step1,
"DeltaSolnTrial",
"dampNewt: Important Entries for "
"Weighted Soln Updates:",
@ -1093,8 +1100,8 @@ namespace Cantera {
*
*/
int NonlinearSolver::solve_nonlinear_problem(int SolnType, double* y_comm,
double* ydot_comm, double CJ,
double time_curr,
double* ydot_comm, doublereal CJ,
doublereal time_curr,
SquareMatrix& jac,
int &num_newt_its,
int &num_linear_solves,
@ -1108,7 +1115,7 @@ namespace Cantera {
bool m_residCurrent = false;
int m = 0;
bool forceNewJac = false;
double s1=1.e30;
doublereal s1=1.e30;
// std::vector<doublereal> y_curr(neq_, 0.0);
std::vector<doublereal> ydot_curr(neq_, 0.0);
@ -1256,7 +1263,7 @@ namespace Cantera {
/*
* Impose max newton iteration
*/
if (num_newt_its > 20) {
if (num_newt_its > maxNewtIts_) {
m = -1;
if (m_print_flag > 1) {
printf("\t\tsolve_nonlinear_problem(): Damped newton unsuccessful (max newts exceeded) sfinal = %g\n", s1);
@ -1355,7 +1362,7 @@ namespace Cantera {
num_linear_solves += m_numTotalLinearSolves;
double time_elapsed = wc.secondsWC();
doublereal time_elapsed = wc.secondsWC();
if (m_print_flag > 1) {
if (m > 0) {
if (NonlinearSolver::m_TurnOffTiming) {
@ -1375,18 +1382,18 @@ namespace Cantera {
*
*/
void NonlinearSolver::
print_solnDelta_norm_contrib(const double * const solnDelta0,
print_solnDelta_norm_contrib(const doublereal * const solnDelta0,
const char * const s0,
const double * const solnDelta1,
const doublereal * const solnDelta1,
const char * const s1,
const char * const title,
const double * const y0,
const double * const y1,
double damp,
const doublereal * const y0,
const doublereal * const y1,
doublereal damp,
int num_entries) {
int i, j, jnum;
bool used;
double dmax0, dmax1, error, rel_norm;
doublereal dmax0, dmax1, error, rel_norm;
printf("\t\t%s currentDamp = %g\n", title, damp);
printf("\t\t I ysolnOld %13s ysolnNewRaw | ysolnNewTrial "
"%10s ysolnNewTrialRaw | solnWeight wtDelSoln wtDelSolnTrial\n", s0, s1);
@ -1442,11 +1449,11 @@ namespace Cantera {
* loss of convergence. Therefore, in practice this routine
* has proved cost-effective.
*/
static inline double subtractRD(double a, double b) {
double diff = a - b;
double d = MIN(fabs(a), fabs(b));
static inline doublereal subtractRD(doublereal a, doublereal b) {
doublereal diff = a - b;
doublereal d = MIN(fabs(a), fabs(b));
d *= 1.0E-14;
double ad = fabs(diff);
doublereal ad = fabs(diff);
if (ad < 1.0E-300) {
diff = 0.0;
}
@ -1467,15 +1474,15 @@ namespace Cantera {
* not have to be computed again.
*
*/
void NonlinearSolver::beuler_jac(SquareMatrix &J, double * const f,
double time_curr, double CJ,
double * const y,
double * const ydot,
void NonlinearSolver::beuler_jac(SquareMatrix &J, doublereal * const f,
doublereal time_curr, doublereal CJ,
doublereal * const y,
doublereal * const ydot,
int num_newt_its)
{
int i, j;
double* col_j;
double ysave, ydotsave, dy;
doublereal ysave, ydotsave, dy;
/*
* Clear the factor flag
*/
@ -1486,7 +1493,7 @@ namespace Cantera {
*/
m_func->evalJacobian(time_curr, delta_t_n, y, ydot, J, f);
#ifdef DEBUG_HKM
//double dddd = J(89, 89);
//doublereal dddd = J(89, 89);
//checkFinite(dddd);
#endif
m_nJacEval++;
@ -1510,7 +1517,7 @@ namespace Cantera {
* deltaY's that are appropriate for calculating the numerical
* derivative.
*/
double *dyVector = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT);
doublereal *dyVector = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT);
m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt));
@ -1546,7 +1553,7 @@ namespace Cantera {
*/
col_j = (double *) J.ptrColumn(j);
col_j = (doublereal *) J.ptrColumn(j);
ysave = y[j];
dy = dyVector[j];
//dy = fmaxx(1.0E-6 * m_ewt[j], fabs(ysave)*1.0E-7);
@ -1565,7 +1572,7 @@ namespace Cantera {
m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp),
JacDelta_ResidEval, j, dy);
m_nfe++;
double diff;
doublereal diff;
for (i = 0; i < neq_; i++) {
diff = subtractRD(m_wksp[i], f[i]);
col_j[i] = diff / dy;
@ -1593,13 +1600,13 @@ namespace Cantera {
* @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr
*/
void NonlinearSolver::
calc_ydot(const int order, const double * const y_curr, double * const ydot_curr)
calc_ydot(const int order, const doublereal * const y_curr, doublereal * const ydot_curr)
{
if (!ydot_curr) {
return;
}
int i;
double c1;
doublereal c1;
switch (order) {
case 0:
case 1: /* First order forward Euler/backward Euler */
@ -1626,7 +1633,7 @@ namespace Cantera {
*
* @return Returns the norm of the value of the amount filtered
*/
double NonlinearSolver::filterNewStep(const double timeCurrent, double * const y_current, double *const ydot_current) {
doublereal NonlinearSolver::filterNewStep(const doublereal timeCurrent, doublereal * const y_current, doublereal *const ydot_current) {
return 0.0;
}
//====================================================================================================================
@ -1661,7 +1668,7 @@ namespace Cantera {
* @param residWts Vector of length neq_
*/
void
NonlinearSolver::getResidWts(double * const residWts) const
NonlinearSolver::getResidWts(doublereal * const residWts) const
{
for (int i = 0; i < neq_; i++) {
residWts[i] = (m_residWts)[i];
@ -1684,7 +1691,7 @@ namespace Cantera {
* 0 Not converged yet
*/
int
NonlinearSolver::convergenceCheck(int dampCode, double s1)
NonlinearSolver::convergenceCheck(int dampCode, doublereal s1)
{
int retn = 0;
if (m_dampBound < 0.9999) {
@ -1743,6 +1750,18 @@ namespace Cantera {
atolk_[i]= atol[i];
}
}
//=====================================================================================================================
// Set the relative tolerances for the solution variables
/*
* Set the relative tolerances used in the calculation
*
* @param rtol single double
*/
void NonlinearSolver::setRtol(const doublereal rtol)
{
rtol_ = rtol;
}
//=====================================================================================================================
}

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

@ -79,7 +79,7 @@ namespace Cantera {
*
* param y vector of the current solution values
*/
void createSolnWeights(const double * const y);
void createSolnWeights(const doublereal * const y);
//! L2 norm of the delta of the solution vector
@ -97,8 +97,8 @@ namespace Cantera {
* @param dampFactor Current value of the damping factor. Defaults to 1.
* only used for printout out a table.
*/
double solnErrorNorm(const double * const delta_y, const char * title = 0, int printLargest = 0,
const double dampFactor = 1.0);
doublereal solnErrorNorm(const doublereal * const delta_y, const char * title = 0, int printLargest = 0,
const doublereal dampFactor = 1.0);
//! L2 norm of the residual of the equation system
/*!
@ -114,8 +114,8 @@ namespace Cantera {
* @param printLargest Number of specific entries to be printed
* @param y Current value of y - only used for printouts
*/
double residErrorNorm(const double * const resid, const char * title = 0, const int printLargest = 0,
const double * const y = 0);
doublereal residErrorNorm(const doublereal * const resid, const char * title = 0, const int printLargest = 0,
const doublereal * const y = 0);
//! Compute the current Residual
/*!
@ -142,8 +142,8 @@ namespace Cantera {
* @param y_curr Current value of the solution vector
* @param ydot_curr Current value of the time derivative of the solution vector
*/
void doResidualCalc(const double time_curr, const int typeCalc, const double * const y_curr,
const double * const ydot_curr);
void doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal * const y_curr,
const doublereal * const ydot_curr);
//! Compute the undamped Newton step
/*!
@ -164,8 +164,8 @@ namespace Cantera {
* @param ydot_current Current value of the solution derivative.
*
*/
void doNewtonSolve(const double time_curr, const double * const y_curr,
const double * const ydot_curr, double* const delta_y,
void doNewtonSolve(const doublereal time_curr, const doublereal * const y_curr,
const doublereal * const ydot_curr, doublereal * const delta_y,
SquareMatrix& jac, int loglevel);
@ -182,7 +182,7 @@ namespace Cantera {
/*!
* @param deltaboundsMagnitudes
*/
void setDeltaBoundsMagnitudes(const double * const deltaBoundsMagnitudes);
void setDeltaBoundsMagnitudes(const doublereal * const deltaBoundsMagnitudes);
//! Bound the step
@ -211,7 +211,7 @@ namespace Cantera {
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*/
double boundStep(const double* const y,
doublereal boundStep(const double* const y,
const double* const step0, const int loglevel);
@ -221,8 +221,8 @@ namespace Cantera {
* @param y_low_bounds Vector of lower bounds
* @param y_high_bounds Vector of high bounds
*/
void setBoundsConstraints(const double * const y_low_bounds,
const double * const y_high_bounds);
void setBoundsConstraints(const doublereal * const y_low_bounds,
const doublereal * const y_high_bounds);
//! Internal function to calculate the time derivative at the new step
@ -231,7 +231,7 @@ namespace Cantera {
* @param y_curr current value of the solution
* @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr
*/
void calc_ydot(const int order, const double * const y_curr, double * const ydot_curr);
void calc_ydot(const int order, const doublereal * const y_curr, doublereal * const ydot_curr);
//! Function called to evaluate the jacobian matrix and the curent
//! residual vector.
@ -239,9 +239,9 @@ namespace Cantera {
*
*
*/
void beuler_jac(SquareMatrix &J, double * const f,
double time_curr, double CJ, double * const y,
double * const ydot, int num_newt_its);
void 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 step
@ -252,7 +252,7 @@ namespace Cantera {
*
* @return Returns the norm of the value of the amount filtered
*/
double filterNewStep(const double timeCurrent, double * const y_current, double * const ydot_current);
doublereal filterNewStep(const doublereal timeCurrent, doublereal * const y_current, doublereal * const ydot_current);
//! Return the factor by which the undamped Newton step 'step0'
@ -272,7 +272,7 @@ namespace Cantera {
*
* @return returns the damping factor
*/
double deltaBoundStep(const double * const y, const double * const step0, const int loglevel);
doublereal deltaBoundStep(const doublereal * const y, const doublereal * const step0, const int loglevel);
//! Find a damping coefficient through a look-ahead mechanism
/*!
@ -292,8 +292,8 @@ namespace Cantera {
* @param step0 Initial step suggested.
* @param y1
*/
int dampStep(const double time_curr, const double* y0,
const double *ydot0, const double* step0,
int dampStep(const doublereal time_curr, const double* y0,
const doublereal *ydot0, const double* step0,
double* const y1, double* const ydot1, double* step1,
double& s1, SquareMatrix& jac,
int& loglevel, bool writetitle,
@ -317,8 +317,8 @@ namespace Cantera {
* -1 Failed convergence
*/
int solve_nonlinear_problem(int SolnType, double* y_comm,
double* ydot_comm, double CJ,
double time_curr,
double* ydot_comm, doublereal CJ,
doublereal time_curr,
SquareMatrix& jac,
int &num_newt_its,
int &num_linear_solves,
@ -333,19 +333,19 @@ namespace Cantera {
/*!
*
*/
void scaleMatrix(SquareMatrix& jac, double* y_comm, double* ydot_comm, double time_curr);
void scaleMatrix(SquareMatrix& jac, double* y_comm, double* ydot_comm, doublereal time_curr);
//! Print solution norm contribution
void
print_solnDelta_norm_contrib(const double * const solnDelta0,
print_solnDelta_norm_contrib(const doublereal * const solnDelta0,
const char * const s0,
const double * const solnDelta1,
const doublereal * const solnDelta1,
const char * const s1,
const char * const title,
const double * const y0,
const double * const y1,
double damp,
const doublereal * const y0,
const doublereal * const y1,
doublereal damp,
int num_entries);
//! Compute the Residual Weights
@ -365,7 +365,7 @@ namespace Cantera {
/*!
* @param residWts Vector of length neq_
*/
void getResidWts(double * const residWts) const;
void getResidWts(doublereal * const residWts) const;
//! Check to see if the nonlinear problem has converged
@ -381,7 +381,7 @@ namespace Cantera {
* The predicted deltaSoln is below 1.0.
* 0 Not converged yet
*/
int convergenceCheck(int dampCode, double s1);
int convergenceCheck(int dampCode, doublereal s1);
//! Set the absolute tolerances for the solution variables
@ -392,6 +392,20 @@ namespace Cantera {
*/
void setAtol(const doublereal * const atol);
//! Set the relative tolerances for the solution variables
/*!
* Set the relative tolerances used in the calculation
*
* @param rtol single double
*/
void setRtol(const doublereal rtol);
//! Set the value of the maximum # of newton iterations
/*!
* @param maxNewtIts Maximum number of newton iterations
*/
void setMaxNewtIts(const int maxNewtIts);
private:
//! Pointer to the residual and jacobian evaluator for the
@ -410,10 +424,6 @@ namespace Cantera {
//! Soln error weights
std::vector<doublereal> m_ewt;
//! Boolean indicating whether a manual delta bounds has been input.
int m_manualDeltaBoundsSet;
@ -462,19 +472,19 @@ namespace Cantera {
std::vector<doublereal> m_residWts;
//! Norm of the residual at the start of each nonlinear iteration
double m_normResid0;
doublereal m_normResid0;
//! Norm of the residual before damping
double m_normResidFRaw;
doublereal m_normResidFRaw;
//! Norm of the solution update created by the iteration in its raw, undamped form.
double m_normSolnFRaw;
doublereal m_normSolnFRaw;
//! Norm of the residual for a trial calculation which may or may not be used
double m_normResidTrial;
doublereal m_normResidTrial;
//! Vector of the norm
double m_normResidPoints[15];
doublereal m_normResidPoints[15];
bool m_resid_scaled;
@ -491,13 +501,13 @@ namespace Cantera {
std::vector<doublereal> m_y_low_bounds;
//! Damping factor imposed by hard bounds and by delta bounds
double m_dampBound;
doublereal m_dampBound;
//! Additional damping factor due to bounds on the residual and solution norms
double m_dampRes;
doublereal m_dampRes;
//! Delta t for the current step
double delta_t_n;
doublereal delta_t_n;
//! Counter for the total number of function evaluations
int m_nfe;
@ -525,6 +535,9 @@ namespace Cantera {
//! Boolean that turns on solution filtering
int filterNewstep;
//! Maximum number of newton iterations
int maxNewtIts_;
//! Jacobian formation method
/*!
* 1 = numerical (default)
@ -540,18 +553,19 @@ namespace Cantera {
* Note, we assume even for steady state problems that the residual
* is a function of a system time.
*/
double time_n;
doublereal time_n;
int m_matrixConditioning;
int m_order;
//! value of the relative tolerance to use in solving the equation set
doublereal rtol_;
//! Base value of the absolute tolerance
doublereal atolBase_;
double * m_ydot_nm1;
doublereal * m_ydot_nm1;
std::vector<doublereal> atolk_;