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

Added in a more robust treatment of handling limits

Browse source
This commit is contained in:
Harry Moffat 2011-01-05 03:55:13 +00:00
parent a36864543a
commit f727cd44af
3 changed files with 52 additions and 9 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -1585,6 +1585,7 @@ namespace Cantera {
{
int i, j;
double* col_j;
int info;
doublereal ysave, ydotsave, dy;
int retn = 1;
/*
@ -1595,9 +1596,12 @@ namespace Cantera {
/********************************************************************
* Call the function to get a jacobian.
*/
m_func->evalJacobian(time_curr, delta_t_n, y, ydot, J, f);
info = m_func->evalJacobian(time_curr, delta_t_n, y, ydot, J, f);
m_nJacEval++;
m_nfe++;
if (info != 1) {
return info;
}
} else {
/*******************************************************************
* Generic algorithm to calculate a numerical Jacobian
@ -1607,8 +1611,11 @@ namespace Cantera {
* current conditions.
*/
m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
m_nfe++;
if (info != 1) {
return info;
}
m_nJacEval++;
/*
@ -1668,9 +1675,14 @@ namespace Cantera {
*/
m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp),
JacDelta_ResidEval, j, dy);
info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp),
JacDelta_ResidEval, j, dy);
m_nfe++;
if (info != 1) {
mdp::mdp_safe_free((void **) &dyVector);
return info;
}
doublereal diff;
for (i = 0; i < neq_; i++) {
diff = subtractRD(m_wksp[i], f[i]);

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

@ -188,7 +188,7 @@ namespace Cantera {
* return:
* 0 Found function
*/
int RootFind::solve(doublereal xmin, doublereal xmax, int itmax, doublereal funcTargetValue, doublereal *xbest) {
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
@ -201,6 +201,8 @@ namespace Cantera {
const char *stre = "RootFind ERROR: ";
const char *strw = "RootFind WARNING: ";
int converged = 0;
int bottomBump = 0;
int topBump = 0;
#ifdef DEBUG_MODE
char fileName[80];
FILE *fp = 0;
@ -236,6 +238,7 @@ namespace Cantera {
#endif
if (xmax <= xmin) {
writelogf("%sxmin and xmax are bad: %g %g\n", stre, xmin, xmax);
funcTargetValue = func(*xbest);
return ROOTFIND_BADINPUT;
}
/*
@ -566,7 +569,18 @@ namespace Cantera {
* Guard against going above xmax or below xmin
*/
if (xnew > xmax) {
xnew = x2 + (xmax - x2) / 2.0;
topBump++;
if (topBump < 5) {
xnew = x2 + (xmax - x2) / 2.0;
} else {
if (x2 == xmax || x1 == xmax) {
// we are here when we are bumping against the top limit.
// No further action is possible
goto done;
} else {
xnew = xmax;
}
}
#ifdef DEBUG_MODE
if (printLvl >= 3) {
fprintf(fp, " | xlimitmax = %-11.5E", xnew);
@ -574,7 +588,18 @@ namespace Cantera {
#endif
}
if (xnew < xmin) {
xnew = x2 + (x2 - xmin) / 2.0;
bottomBump++;
if (bottomBump < 5) {
xnew = x2 + (x2 - xmin) / 2.0;
} else {
if (x2 == xmin || x1 == xmin) {
// we are here when we are bumping against the bottom limit.
// No further action is possible
goto done;
} else {
xnew = xmin;
}
}
#ifdef DEBUG_MODE
if (printLvl >= 3) {
fprintf(fp, " | xlimitmin = %-11.5E", xnew);
@ -658,6 +683,8 @@ namespace Cantera {
}
its++;
} while (! converged && its < itmax);
done:
if (converged) {
if (printLvl >= 1) {
writelogf("RootFind success: convergence achieved\n");
@ -679,11 +706,13 @@ namespace Cantera {
#endif
}
*xbest = x2;
funcTargetValue = f2 + m_funcTargetValue;
#ifdef DEBUG_MODE
if (printLvl >= 3) {
fclose(fp);
}
#endif
return retn;
}
//================================================================================================

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

@ -76,7 +76,9 @@ namespace Cantera {
* @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 funcTargetValue
* Value of \f$ f_o \f$ in the equation.
* On return, it contains the value of the function actually obtained.
* @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.
@ -86,7 +88,7 @@ namespace Cantera {
* -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);
int solve(doublereal xmin, doublereal xmax, int itmax, doublereal &funcTargetValue, doublereal *xbest);
//! Return the function value