Revamped the rootfinder algorithm.
Behavior is substantially different. A summary table is now written out at the end of the solve for some print levels.
This commit is contained in:
Harry Moffat
parent
445adce8fc
commit
d251da7a84
2 changed files with 510 additions and 136 deletions
|
|
@ -18,7 +18,7 @@
|
|||
#ifdef DEBUG_MODE
|
||||
#include "mdp_allo.h"
|
||||
#endif
|
||||
|
||||
#include "stringUtils.h"
|
||||
/* Standard include files */
|
||||
|
||||
#include <cstdio>
|
||||
|
|
@ -48,6 +48,13 @@ namespace Cantera {
|
|||
#define DSIGN(x) (( (x) == (0.0) ) ? (0.0) : ( ((x) > 0.0) ? 1.0 : -1.0 ))
|
||||
#endif
|
||||
|
||||
#ifdef SWAP
|
||||
#undef SWAP
|
||||
#endif
|
||||
#ifndef SWAP
|
||||
#define SWAP(x1, x2, tmp) ((tmp) = (x2), (x2) = (x1), (x1) = (tmp))
|
||||
#endif
|
||||
|
||||
// turn on debugging for now
|
||||
#ifndef DEBUG_MODE
|
||||
#define DEBUG_MODE
|
||||
|
|
@ -275,17 +282,20 @@ namespace Cantera {
|
|||
// 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
|
||||
* are close to one another. This is defined as being within factor * delXMeaningful() of each other.
|
||||
*
|
||||
* @param x2 First number
|
||||
* @param x2 second 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) const {
|
||||
bool RootFind::theSame(doublereal x2, doublereal x1, doublereal factor) const {
|
||||
doublereal x = fabs(x2) + fabs(x1);
|
||||
doublereal deltaX = delXMeaningful(x);
|
||||
if (fabs(x2 - x1) < deltaX) {
|
||||
doublereal deltaXSmall = factor * deltaX;
|
||||
deltaXSmall = MAX(deltaXSmall , x * 1.0E-15);
|
||||
if (fabs(x2 - x1) < deltaXSmall) {
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
|
|
@ -322,9 +332,10 @@ namespace Cantera {
|
|||
char fileName[80];
|
||||
FILE *fp = 0;
|
||||
#endif
|
||||
int doFinalFuncCall = 0;
|
||||
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 = ROOTFIND_FAILEDCONVERGENCE;
|
||||
int foundPosF = 0;
|
||||
|
|
@ -335,11 +346,15 @@ namespace Cantera {
|
|||
doublereal xNegF = 0.0;
|
||||
doublereal fNegF = -1.0E300;
|
||||
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 x0 = 0.0, f0 = 0.0, xDelMin;
|
||||
doublereal sgn;
|
||||
|
||||
|
||||
doublereal dtmp;
|
||||
doublereal fnoise = 0.0;
|
||||
rfHistory_.clear();
|
||||
rfTable rfT;
|
||||
rfT.clear();
|
||||
rfT.reasoning = "First Point: ";
|
||||
|
||||
callNum++;
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
|
|
@ -397,10 +412,14 @@ namespace Cantera {
|
|||
x1 = *xbest;
|
||||
if (x1 < xmin || x1 > xmax) {
|
||||
x1 = (xmin + xmax) / 2.0;
|
||||
rfT.reasoning += " x1 set middle between xmin and xmax because entrance is outside bounds.";
|
||||
} else {
|
||||
rfT.reasoning += " x1 set to entrance x.";
|
||||
}
|
||||
|
||||
x_maxTried_ = x1;
|
||||
x_minTried_ = x1;
|
||||
int its = 1;
|
||||
f1 = func(x1);
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
|
|
@ -410,45 +429,62 @@ namespace Cantera {
|
|||
}
|
||||
#endif
|
||||
|
||||
|
||||
if (f1 == 0.0) {
|
||||
*xbest = x1;
|
||||
return 0;
|
||||
} else if (f1 > 0.0) {
|
||||
} else if (f1 > fnoise) {
|
||||
foundPosF = 1;
|
||||
xPosF = x1;
|
||||
fPosF = f1;
|
||||
} else {
|
||||
} else if (f1 < -fnoise) {
|
||||
foundNegF = 1;
|
||||
xNegF = x1;
|
||||
fNegF = x1;
|
||||
fNegF = f1;
|
||||
}
|
||||
|
||||
rfT.its = its;
|
||||
rfT.TP_its = 0;
|
||||
rfT.xval = x1;
|
||||
rfT.fval = f1;
|
||||
rfT.foundPos = foundPosF;
|
||||
rfT.foundNeg = foundNegF;
|
||||
rfT.deltaXConverged = m_rtolx * (fabs(x1) + 0.001);
|
||||
rfT.deltaFConverged = fabs(f1) * m_rtolf;
|
||||
rfT.delX = xmax - xmin;
|
||||
rfHistory_.push_back(rfT);
|
||||
rfT.clear();
|
||||
|
||||
/*
|
||||
* Now, this is actually a tricky part of the algorithm - Find the x value for
|
||||
* the second point. It's tricky because we don't have a valid idea of the scale of x yet
|
||||
*
|
||||
*/
|
||||
rfT.reasoning = "Second Point: ";
|
||||
if (x1 == 0.0) {
|
||||
x2 = x1 + 0.01 * DeltaXnorm_;
|
||||
rfT.reasoning += "Set by DeltaXnorm_";
|
||||
} else {
|
||||
x2 = x1 * 1.0001;
|
||||
rfT.reasoning += "Set slightly higher.";
|
||||
}
|
||||
if (x2 > xmax) {
|
||||
x2 = x1 - 0.01 * DeltaXnorm_;
|
||||
rfT.reasoning += " - But adjusted to be within bounds";
|
||||
}
|
||||
|
||||
/*
|
||||
* Find the second function value f2 = func(x2), Process it
|
||||
*/
|
||||
deltaX2 = x2 - x1;
|
||||
f2 = func(x2);
|
||||
its++;
|
||||
f2 = func(x2);
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
print_funcEval(fp, x2, f2, its);
|
||||
fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", -1, 0, x2, f2);
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
/*
|
||||
* Calculate the norm of the function, this is the nominal value of f. We try
|
||||
* to reduce the nominal value of f by rtolf, this is the main convergence requirement.
|
||||
|
|
@ -458,23 +494,33 @@ namespace Cantera {
|
|||
} else {
|
||||
fnorm = 0.5*(fabs(f1) + fabs(f2)) + fabs(m_funcTargetValue) + m_atolf;
|
||||
}
|
||||
fnoise = 1.0E-100;
|
||||
|
||||
if (f2 == 0.0) {
|
||||
*xbest = x2;
|
||||
return ROOTFIND_SUCCESS;
|
||||
} else if (f2 > 0.0) {
|
||||
|
||||
if (f2 > fnoise) {
|
||||
if (!foundPosF) {
|
||||
foundPosF = 1;
|
||||
xPosF = x2;
|
||||
fPosF = x2;
|
||||
fPosF = f2;
|
||||
}
|
||||
} else {
|
||||
} else if (f2 < - fnoise) {
|
||||
if (!foundNegF) {
|
||||
foundNegF = 1;
|
||||
xNegF = x2;
|
||||
fNegF = f2;
|
||||
}
|
||||
} else if (f2 == 0.0) {
|
||||
*xbest = x2;
|
||||
return ROOTFIND_SUCCESS;
|
||||
}
|
||||
rfT.its = its;
|
||||
rfT.TP_its = 0;
|
||||
rfT.xval = x2;
|
||||
rfT.fval = f2;
|
||||
rfT.foundPos = foundPosF;
|
||||
rfT.foundNeg = foundNegF;
|
||||
|
||||
|
||||
/*
|
||||
* See if we have already achieved a straddle
|
||||
*/
|
||||
|
|
@ -483,7 +529,7 @@ namespace Cantera {
|
|||
if (xPosF > xNegF) posStraddle = 1;
|
||||
else posStraddle = 0;
|
||||
}
|
||||
bool doQuad = false;
|
||||
|
||||
bool useNextStrat = false;
|
||||
bool slopePointingToHigher = true;
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
|
|
@ -499,8 +545,13 @@ namespace Cantera {
|
|||
printf(" RootFind: we are here x2 = %g x1 = %g\n", x2, x1);
|
||||
}
|
||||
#endif
|
||||
|
||||
doublereal delXtmp = deltaXControlled(x2, x1);
|
||||
slope = (f2 - f1) / delXtmp;
|
||||
rfT.slope = slope;
|
||||
rfHistory_.push_back(rfT);
|
||||
rfT.clear();
|
||||
rfT.reasoning = "";
|
||||
if (fabs(slope) <= 1.0E-100) {
|
||||
if (printLvl >= 2) {
|
||||
writelogf("%s functions evals produced the same result, %g, at %g and %g\n",
|
||||
|
|
@ -509,6 +560,7 @@ namespace Cantera {
|
|||
xnew = x2 + DeltaXnorm_;
|
||||
slopePointingToHigher = true;
|
||||
useNextStrat = true;
|
||||
rfT.reasoning += "Slope is close to zero. ";
|
||||
} else {
|
||||
useNextStrat = false;
|
||||
xnew = x2 - f2 / slope;
|
||||
|
|
@ -517,6 +569,7 @@ namespace Cantera {
|
|||
} else {
|
||||
slopePointingToHigher = false;
|
||||
}
|
||||
rfT.reasoning += "Slope is good. ";
|
||||
}
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
|
|
@ -530,13 +583,18 @@ namespace Cantera {
|
|||
if (!foundStraddle) {
|
||||
if (fabs(xnew - x2) > DeltaXMax_) {
|
||||
useNextStrat = true;
|
||||
rfT.reasoning += "Too large change in xnew from slope. ";
|
||||
}
|
||||
if (fabs(deltaXnew) < fabs(deltaX2)) {
|
||||
deltaXnew = 1.2 * deltaXnew;
|
||||
xnew = x2 + deltaXnew;
|
||||
}
|
||||
}
|
||||
/*
|
||||
* If the slope can't be trusted using a different strategy for picking the next point
|
||||
*/
|
||||
if (useNextStrat) {
|
||||
rfT.reasoning += "Using DeltaXnorm, " + fp2str(DeltaXnorm_) + " and FuncIsGenerallyIncreasing hints. ";
|
||||
if (f2 < 0.0) {
|
||||
if (FuncIsGenerallyIncreasing_) {
|
||||
if (slopePointingToHigher) {
|
||||
|
|
@ -580,53 +638,24 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Do a quadratic fit -> Note this algorithm seems
|
||||
* to work OK. The quadratic approximation doesn't kick in until
|
||||
* the end of the run, when it becomes reliable.
|
||||
* Here, if we have a straddle, we purposefully overshoot the smaller side by 5%. Yes it does lead to
|
||||
* more iterations. However, we're interested in bounding x, and not just doing Newton's method.
|
||||
*/
|
||||
if (its > 0 && doQuad) {
|
||||
c[0] = 1.; c[1] = 1.; c[2] = 1.;
|
||||
c[3] = x0; c[4] = x1; c[5] = x2;
|
||||
c[6] = SQUARE(x0); c[7] = SQUARE(x1); c[8] = SQUARE(x2);
|
||||
f[0] = - f0; f[1] = - f1; f[2] = - f2;
|
||||
int rrr = smlequ(c, 3, 3, f, 1);
|
||||
if (rrr == 1) goto QUAD_BAIL;
|
||||
root = f[1]* f[1] - 4.0 * f[0] * f[2];
|
||||
if (root >= 0.0) {
|
||||
xn1 = (- f[1] + sqrt(root)) / (2.0 * f[2]);
|
||||
xn2 = (- f[1] - sqrt(root)) / (2.0 * f[2]);
|
||||
if (fabs(xn2 - x2) < fabs(xn1 - x2) && xn2 > 0.0 ) xquad = xn2;
|
||||
else xquad = xn1;
|
||||
theta = fabs(xquad - xnew) / fabs(xnew - x2);
|
||||
theta = MIN(1.0, theta);
|
||||
xnew = theta * xnew + (1.0 - theta) * xquad;
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
if (theta != 1.0) {
|
||||
fprintf(fp, " | xquad = %-11.5E", xnew);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
} else {
|
||||
/*
|
||||
* Pick out situations where the convergence may be
|
||||
* accelerated.
|
||||
*/
|
||||
if ((DSIGN(xnew - x2) == DSIGN(x2 - x1)) &&
|
||||
(DSIGN(x2 - x1) == DSIGN(x1 - x0)) ) {
|
||||
xnew += xnew - x2;
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
fprintf(fp, " | xquada = %-11.5E", xnew);
|
||||
}
|
||||
#endif
|
||||
if (foundStraddle) {
|
||||
double delta = fabs(x2 - x1);
|
||||
if (fabs(xnew - x1) < .01 * delta) {
|
||||
xnew = x1 + 0.01 * (x2 - x1);
|
||||
} else if (fabs(xnew - x2) < .01 * delta) {
|
||||
xnew = x1 + 0.01 * (x2 - x1);
|
||||
} else if ((xnew > x1 && xnew < x2) || (xnew < x1 && xnew > x2)) {
|
||||
if (fabs(xnew - x1) < fabs(x2 - xnew)) {
|
||||
xnew = x1 + 20./19. * (xnew - x1);
|
||||
} else {
|
||||
xnew = x2 + 20./19. * (xnew - x2);
|
||||
}
|
||||
}
|
||||
}
|
||||
QUAD_BAIL: ;
|
||||
|
||||
/*
|
||||
* OK, we have an estimate xnew.
|
||||
*
|
||||
|
|
@ -637,7 +666,7 @@ namespace Cantera {
|
|||
/*
|
||||
* If we are doing a jump in between the two previous points, make sure
|
||||
* the new trial is no closer that 10% of the distances between x2-x1 to
|
||||
* any of the original points.
|
||||
* any of the original points. This is an important part of finding a good bound.
|
||||
*/
|
||||
xDelMin = fabs(x2 - x1) / 10.;
|
||||
if (fabs(xnew - x1) < xDelMin) {
|
||||
|
|
@ -674,13 +703,37 @@ namespace Cantera {
|
|||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
fprintf(fp, " | xlimitsize = %-11.5E", xnew);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
/*
|
||||
* If we are doing a jump outside the two previous points, make sure
|
||||
* the new trial is no closer that 10% of the distances between x2-x1 to
|
||||
* any of the original points. This is an important part of finding a good bound.
|
||||
*/
|
||||
xDelMin = 0.1 * fabs(x2 - x1);
|
||||
if (fabs(xnew - x2) < xDelMin) {
|
||||
xnew = x2 + DSIGN(xnew - x2) * xDelMin;
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
fprintf(fp, " | x10%% = %-11.5E", xnew);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
if (fabs(xnew - x1) < xDelMin) {
|
||||
xnew = x1 + DSIGN(xnew - x1) * xDelMin;
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
fprintf(fp, " | x10%% = %-11.5E", xnew);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* HKM -> Not sure this section is needed
|
||||
*/
|
||||
if (foundStraddle) {
|
||||
#ifdef DEBUG_MODE
|
||||
slope = xnew;
|
||||
double xorig = xnew;
|
||||
#endif
|
||||
if (posStraddle) {
|
||||
if (f2 > 0.0) {
|
||||
|
|
@ -717,7 +770,7 @@ namespace Cantera {
|
|||
}
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
if (slope != xnew) {
|
||||
if (xorig != xnew) {
|
||||
fprintf(fp, " | xstraddle = %-11.5E", xnew);
|
||||
}
|
||||
}
|
||||
|
|
@ -735,6 +788,8 @@ namespace Cantera {
|
|||
sgn = -1.0;
|
||||
}
|
||||
deltaXnew = 1.2 * delXMeaningful(xnew) * sgn;
|
||||
rfT.reasoning += "Enforcing minimum stepsize from " + fp2str(xnew - x2) +
|
||||
" to " + fp2str(deltaXnew);
|
||||
xnew = x2 + deltaXnew;
|
||||
}
|
||||
}
|
||||
|
|
@ -746,20 +801,18 @@ namespace Cantera {
|
|||
topBump++;
|
||||
if (topBump < 3) {
|
||||
xnew = x2 + (xmax - x2) / 2.0;
|
||||
rfT.reasoning += ("xval reduced to " + fp2str(xnew) + " because predicted xnew was above max value of " + fp2str(xmax));
|
||||
} else {
|
||||
if (x2 == xmax || x1 == xmax) {
|
||||
// we are here when we are bumping against the top limit.
|
||||
// No further action is possible
|
||||
if (xnew > xmax) {
|
||||
slope = (f2 - f1) / delXtmp;
|
||||
xnew = x2 - f2 / slope;
|
||||
if (xnew > xmax) {
|
||||
retn = ROOTFIND_SOLNHIGHERTHANXMAX;
|
||||
*xbest = xnew;
|
||||
}
|
||||
}
|
||||
retn = ROOTFIND_SOLNHIGHERTHANXMAX;
|
||||
*xbest = xnew;
|
||||
rfT.slope = slope;
|
||||
rfT.reasoning += "Giving up because we're at xmax and xnew point higher: " + fp2str(xnew);
|
||||
goto done;
|
||||
} else {
|
||||
rfT.reasoning += "xval reduced from " + fp2str(xnew) + " to the max value, " + fp2str(xmax);
|
||||
xnew = xmax;
|
||||
}
|
||||
}
|
||||
|
|
@ -772,21 +825,20 @@ namespace Cantera {
|
|||
if (xnew < xmin) {
|
||||
bottomBump++;
|
||||
if (bottomBump < 3) {
|
||||
rfT.reasoning += ("xnew increased from " + fp2str(xnew) +" to " + fp2str(x2 - (x2 - xmin) / 2.0) +
|
||||
" because above min value of " + fp2str(xmin));
|
||||
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
|
||||
if (xnew < xmin) {
|
||||
slope = (f2 - f1) / delXtmp;
|
||||
xnew = x2 - f2 / slope;
|
||||
if (xnew < xmin) {
|
||||
retn = ROOTFIND_SOLNLOWERTHANXMIN;
|
||||
*xbest = xnew;
|
||||
}
|
||||
}
|
||||
retn = ROOTFIND_SOLNLOWERTHANXMIN;
|
||||
*xbest = xnew;
|
||||
rfT.slope = slope;
|
||||
rfT.reasoning = "Giving up because we're already at xmin and xnew points lower: " + fp2str(xnew);
|
||||
goto done;
|
||||
} else {
|
||||
} else {
|
||||
rfT.reasoning += "xval increased from " + fp2str(xnew) + " to the min value, " + fp2str(xmin);
|
||||
xnew = xmin;
|
||||
}
|
||||
}
|
||||
|
|
@ -796,10 +848,10 @@ namespace Cantera {
|
|||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
its++;
|
||||
fnew = func(xnew);
|
||||
CRdenom = MAX(fabs(fnew), MAX(fabs(f2), MAX(fabs(f1), fnorm)));
|
||||
CRnew = sqrt(fabs(fnew) / CRdenom);
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
fprintf(fp,"\n");
|
||||
|
|
@ -807,7 +859,9 @@ namespace Cantera {
|
|||
fprintf(fp, "%-5d %-5d %-15.5E %-15.5E", its, 0, xnew, fnew);
|
||||
}
|
||||
#endif
|
||||
|
||||
rfT.xval = xnew;
|
||||
rfT.fval = fnew;
|
||||
rfT.its = its;
|
||||
if (foundStraddle) {
|
||||
if (posStraddle) {
|
||||
if (fnew > 0.0) {
|
||||
|
|
@ -837,18 +891,20 @@ namespace Cantera {
|
|||
}
|
||||
|
||||
if (! foundStraddle) {
|
||||
if (fnew > 0.0) {
|
||||
if (fnew > fnoise) {
|
||||
if (!foundPosF) {
|
||||
foundPosF = 1;
|
||||
rfT.foundPos = 1;
|
||||
xPosF = xnew;
|
||||
fPosF = fnew;
|
||||
foundStraddle = 1;
|
||||
if (xPosF > xNegF) posStraddle = 1;
|
||||
else posStraddle = 0;
|
||||
}
|
||||
} else {
|
||||
} else if (fnew < - fnoise) {
|
||||
if (!foundNegF) {
|
||||
foundNegF = 1;
|
||||
rfT.foundNeg = 1;
|
||||
xNegF = xnew;
|
||||
fNegF = fnew;
|
||||
foundStraddle = 1;
|
||||
|
|
@ -860,15 +916,14 @@ namespace Cantera {
|
|||
|
||||
x0 = x1;
|
||||
f0 = f1;
|
||||
CR0 = CR1;
|
||||
x1 = x2;
|
||||
f1 = f2;
|
||||
CR1 = CR2;
|
||||
|
||||
x2 = xnew;
|
||||
f2 = fnew;
|
||||
|
||||
/*
|
||||
* As we go on to new data points, we make sure that
|
||||
* As we go on to new data points, we make sure that
|
||||
* we have the best straddle of the solution with the choice of F1 and F2 when
|
||||
* we do have a straddle to work with.
|
||||
*/
|
||||
|
|
@ -947,64 +1002,190 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
AssertThrow((f1 * f2 <= 0.0), "F1 and F2 aren't bounding");
|
||||
}
|
||||
|
||||
deltaX1 = deltaX2;
|
||||
deltaX2 = deltaXnew;
|
||||
CR2 = CRnew;
|
||||
deltaXnew = x2 - x1;
|
||||
deltaXConverged_ = 0.5 * deltaXConverged_ + 0.5 * (m_rtolx * 0.5 * (fabs(x2) + fabs(x1)) + m_atolx);
|
||||
if (fabs(fnew / fnorm) < m_rtolf) {
|
||||
if (deltaX2 < deltaXConverged_ && deltaXnew < deltaXConverged_) {
|
||||
converged = 1;
|
||||
rfT.deltaXConverged = deltaXConverged_;
|
||||
rfT.deltaFConverged = fnorm * m_rtolf;
|
||||
if (foundStraddle) {
|
||||
rfT.delX = MAX(fabs(deltaX2), fabs(deltaXnew));
|
||||
} else {
|
||||
rfT.delX = MAX(fabs(deltaX2), fabs(deltaXnew));
|
||||
if (x2 < x1) {
|
||||
rfT.delX = MAX(rfT.delX, x2 - xmin);
|
||||
} else {
|
||||
rfT.delX = MAX(rfT.delX, xmax - x2);
|
||||
}
|
||||
if (fabs(slope) > 1.0E-100) {
|
||||
double xdels = fabs(fnew / slope);
|
||||
if (xdels < deltaXConverged_ * 0.5) {
|
||||
converged = 1;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
/*
|
||||
* Check for excess convergence in the x coordinate
|
||||
* Section To Determine CONVERGENCE criteria
|
||||
*/
|
||||
if (foundStraddle) {
|
||||
doublereal denom = fabs(x1) + fabs(x2);
|
||||
if (denom < 1.0E-200) {
|
||||
retn = ROOTFIND_FAILEDCONVERGENCE;
|
||||
converged = true;
|
||||
doFinalFuncCall = 0;
|
||||
if ((fabs(fnew / fnorm) < m_rtolf) && foundStraddle) {
|
||||
if (fabs(deltaX2) < deltaXConverged_ && fabs(deltaXnew) < deltaXConverged_) {
|
||||
converged = 1;
|
||||
rfT.reasoning += "NormalConvergence";
|
||||
retn = ROOTFIND_SUCCESS;
|
||||
}
|
||||
|
||||
else if (fabs(slope) > 1.0E-100) {
|
||||
double xdels = fabs(fnew / slope);
|
||||
if (xdels < deltaXConverged_ * 0.3) {
|
||||
converged = 1;
|
||||
rfT.reasoning += "NormalConvergence-SlopelimitsDelX";
|
||||
doFinalFuncCall = 1;
|
||||
retn = ROOTFIND_SUCCESS;
|
||||
}
|
||||
}
|
||||
if (theSame(x2, x1)) {
|
||||
converged = true;
|
||||
|
||||
|
||||
/*
|
||||
* Check for excess convergence in the x coordinate
|
||||
*/
|
||||
if (!converged) {
|
||||
if (foundStraddle) {
|
||||
doublereal denom = fabs(x1 - x2);
|
||||
if (denom < 1.0E-200) {
|
||||
retn = ROOTFIND_FAILEDCONVERGENCE;
|
||||
converged = true;
|
||||
rfT.reasoning += "ConvergenceFZero but X1X2Identical";
|
||||
}
|
||||
if (theSame(x2, x1, 1.0E-2)) {
|
||||
converged = true;
|
||||
rfT.reasoning += " ConvergenceF and XSame";
|
||||
retn = ROOTFIND_SUCCESS;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
/*
|
||||
* We are here when F is not converged, but we may want to end anyway
|
||||
*/
|
||||
if (!converged) {
|
||||
if (foundStraddle) {
|
||||
doublereal denom = fabs(x1 - x2);
|
||||
if (denom < 1.0E-200) {
|
||||
retn = ROOTFIND_FAILEDCONVERGENCE;
|
||||
converged = true;
|
||||
rfT.reasoning += "FNotConverged but X1X2Identical";
|
||||
}
|
||||
/*
|
||||
* The premise here is that if x1 and x2 get close to one another,
|
||||
* then the accuracy of the calculation gets destroyed.
|
||||
*/
|
||||
if (theSame(x2, x1, 1.0E-5)) {
|
||||
converged = true;
|
||||
retn = ROOTFIND_SUCCESS_XCONVERGENCEONLY;
|
||||
rfT.reasoning += "FNotConverged but XSame";
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
its++;
|
||||
} while (! converged && its < itmax);
|
||||
|
||||
done:
|
||||
if (converged) {
|
||||
retn = ROOTFIND_SUCCESS;
|
||||
if (fabs(f1) < 2.0 * fabs(f2)) {
|
||||
slope = (f2 - f1) / (x2 - x1);
|
||||
xnew = x2 - f2 / slope;
|
||||
rfT.slope = slope;
|
||||
rfHistory_.push_back(rfT);
|
||||
rfT.clear();
|
||||
rfT.its = its;
|
||||
AssertThrow((f1 * f2 <= 0.0), "F1 and F2 aren't bounding");
|
||||
|
||||
double x_fpos = x2;
|
||||
double x_fneg = x1;
|
||||
if (f2 < 0.0) {
|
||||
x_fpos = x1;
|
||||
x_fneg = x2;
|
||||
}
|
||||
rfT.delX = fabs(x_fpos - x_fneg);
|
||||
if (doFinalFuncCall || (fabs(f1) < 2.0 * fabs(f2))) {
|
||||
double delXtmp = deltaXControlled(x2, x1);
|
||||
slope = (f2 - f1) / delXtmp;
|
||||
xnew = x2 - f2 / slope;
|
||||
its++;
|
||||
fnew = func(xnew);
|
||||
if (fabs(fnew) < fabs(f2)) {
|
||||
if (fnew > 0.0) {
|
||||
if (fabs(xnew - x_fneg) < fabs(x_fpos - x_fneg)) {
|
||||
x_fpos = xnew;
|
||||
rfT.delX = fabs(xnew - x_fneg);
|
||||
}
|
||||
} else {
|
||||
if (fabs(xnew - x_fpos) < fabs(x_fpos - x_fneg)) {
|
||||
x_fneg = xnew;
|
||||
rfT.delX = fabs(xnew - x_fpos);
|
||||
}
|
||||
}
|
||||
rfT.its = its;
|
||||
if (fabs(fnew) < fabs(f2) && (fabs(fnew) < fabs(f1))) {
|
||||
*xbest = xnew;
|
||||
if (doFinalFuncCall) {
|
||||
rfT.reasoning += "CONVERGENCE: Another Evaluation Requested";
|
||||
rfT.delX = fabs(xnew - x2);
|
||||
} else {
|
||||
rfT.reasoning += "CONVERGENCE: Another Evaluation done because f1 < f2";
|
||||
rfT.delX = fabs(xnew - x1);
|
||||
}
|
||||
rfT.fval = fnew;
|
||||
rfT.xval = xnew;
|
||||
x2 = xnew;
|
||||
f2 = fnew;
|
||||
} else if (fabs(f1) < fabs(f2)) {
|
||||
rfT.its = its;
|
||||
rfT.xval = xnew;
|
||||
rfT.fval = fnew;
|
||||
|
||||
rfT.slope = slope;
|
||||
rfT.reasoning += "CONVERGENCE: Another Evaluation not as good as Second Point ";
|
||||
rfHistory_.push_back(rfT);
|
||||
rfT.clear();
|
||||
rfT.its = its;
|
||||
SWAP(f1, f2, dtmp);
|
||||
SWAP(x1, x2, dtmp);
|
||||
*xbest = x2;
|
||||
if (fabs(fnew) < fabs(f1)) {
|
||||
if (f1 * fnew > 0.0) {
|
||||
SWAP(f1, fnew, dtmp);
|
||||
SWAP(x1, xnew, dtmp);
|
||||
}
|
||||
}
|
||||
|
||||
rfT.its = its;
|
||||
rfT.xval = *xbest;
|
||||
rfT.fval = f2;
|
||||
rfT.delX = fabs(x_fpos - x_fneg);
|
||||
rfT.reasoning += "CONVERGENCE: NormalEnding -> Second point used";
|
||||
} else {
|
||||
rfT.its = its;
|
||||
rfT.xval = xnew;
|
||||
rfT.fval = fnew;
|
||||
|
||||
rfT.slope = slope;
|
||||
rfT.reasoning += "CONVERGENCE: Another Evaluation not as good as First Point ";
|
||||
rfHistory_.push_back(rfT);
|
||||
rfT.clear();
|
||||
rfT.its = its;
|
||||
*xbest = x2;
|
||||
rfT.xval = *xbest;
|
||||
rfT.fval = f2;
|
||||
rfT.delX = fabs(x_fpos - x_fneg);
|
||||
rfT.reasoning += "CONVERGENCE: NormalEnding -> Last point used";
|
||||
}
|
||||
if (fabs(f1) < fabs(f2)) {
|
||||
x2 = x1;
|
||||
f2 = f1;
|
||||
*xbest = x1;
|
||||
fnew = func(x2);
|
||||
}
|
||||
} else {
|
||||
|
||||
*xbest = x2;
|
||||
|
||||
rfT.xval = *xbest;
|
||||
rfT.fval = f2;
|
||||
rfT.delX = fabs(x2 - x1);
|
||||
rfT.reasoning += "CONVERGENCE: NormalEnding -> Last point used";
|
||||
}
|
||||
|
||||
|
||||
|
||||
funcTargetValue = f2 + m_funcTargetValue;
|
||||
rfT.slope = slope;
|
||||
|
||||
if (printLvl >= 1) {
|
||||
writelogf("RootFind success: convergence achieved\n");
|
||||
|
|
@ -1013,35 +1194,50 @@ namespace Cantera {
|
|||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
fprintf(fp, " | RootFind success in %d its, fnorm = %g\n", its, fnorm);
|
||||
}
|
||||
#endif
|
||||
#endif
|
||||
rfHistory_.push_back(rfT);
|
||||
} else {
|
||||
rfT.reasoning = "FAILED CONVERGENCE ";
|
||||
rfT.slope = slope;
|
||||
rfT.its = its;
|
||||
if (retn == ROOTFIND_SOLNHIGHERTHANXMAX) {
|
||||
if (printLvl >= 1) {
|
||||
writelogf("RootFind ERROR: Soln probably lies higher than xmax, %g: best guess = %g\n", xmax, *xbest);
|
||||
}
|
||||
rfT.reasoning += "Soln probably lies higher than xmax, " + fp2str(xmax) + ": best guess = " + fp2str(*xbest);
|
||||
} else if (retn == ROOTFIND_SOLNLOWERTHANXMIN) {
|
||||
if (printLvl >= 1) {
|
||||
writelogf("RootFind ERROR: Soln probably lies lower than xmin, %g: best guess = %g\n", xmin, *xbest);
|
||||
}
|
||||
rfT.reasoning += "Soln probably lies lower than xmin, " + fp2str(xmin) + ": best guess = " + fp2str(*xbest);
|
||||
} else {
|
||||
retn = ROOTFIND_FAILEDCONVERGENCE;
|
||||
if (printLvl >= 1) {
|
||||
writelogf("RootFind ERROR: maximum iterations exceeded without convergence, cause unknown\n");
|
||||
}
|
||||
rfT.reasoning += "Maximum iterations exceeded without convergence, cause unknown";
|
||||
}
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
fprintf(fp, "\nRootFind failure in %d its\n", its);
|
||||
}
|
||||
#endif
|
||||
|
||||
*xbest = x2;
|
||||
funcTargetValue = f2 + m_funcTargetValue;
|
||||
rfT.xval = *xbest;
|
||||
rfT.fval = f2;
|
||||
rfHistory_.push_back(rfT);
|
||||
}
|
||||
*xbest = x2;
|
||||
funcTargetValue = f2 + m_funcTargetValue;
|
||||
#ifdef DEBUG_MODE
|
||||
if (printLvl >= 3 && writeLogAllowed_) {
|
||||
fclose(fp);
|
||||
}
|
||||
#endif
|
||||
|
||||
if (printLvl >= 2) {
|
||||
printTable();
|
||||
}
|
||||
|
||||
return retn;
|
||||
}
|
||||
|
|
@ -1173,4 +1369,22 @@ namespace Cantera {
|
|||
specifiedDeltaXMax_ = 1;
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
||||
//====================================================================================================================
|
||||
void RootFind::printTable() {
|
||||
printf("\t----------------------------------------------------------------------------------------------------------------------------------------\n");
|
||||
printf("\t RootFinder Summary table: \n");
|
||||
printf("\t FTarget = %g\n", m_funcTargetValue);
|
||||
printf("\t Iter | xval delX deltaXConv | slope | foundP foundN| F - F_targ deltaFConv | Reasoning\n");
|
||||
printf("\t----------------------------------------------------------------------------------------------------------------------------------------\n");
|
||||
for (int i = 0; i < (int) rfHistory_.size(); i++) {
|
||||
struct rfTable rfT = rfHistory_[i];
|
||||
printf("\t %3d |%- 17.11E %- 13.7E %- 13.7E |%- 13.5E| %3d %3d | %- 12.5E %- 12.5E | %s \n",
|
||||
rfT.its, rfT.xval, rfT.delX, rfT.deltaXConverged, rfT.slope, rfT.foundPos, rfT.foundNeg, rfT.fval,
|
||||
rfT.deltaFConverged, (rfT.reasoning).c_str());
|
||||
}
|
||||
printf("\t----------------------------------------------------------------------------------------------------------------------------------------\n");
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -29,6 +29,14 @@ namespace Cantera {
|
|||
|
||||
//! This means that the root solver was a success
|
||||
#define ROOTFIND_SUCCESS 0
|
||||
//! This return value means that the root finder resolved a solution in the x coordinate
|
||||
//! However, convergence in F was not achieved.
|
||||
/*!
|
||||
* A common situation for this to happen is that f(x) is discontinuous about f(x) = f_0,
|
||||
* where we seek the x where the function is equal to f_0. f(x) spans the
|
||||
* f_0 while not being equal to f_0 anywhere.
|
||||
*/
|
||||
#define ROOTFIND_SUCCESS_XCONVERGENCEONLY 1
|
||||
//! 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
|
||||
|
|
@ -40,12 +48,96 @@ namespace Cantera {
|
|||
*
|
||||
*/
|
||||
#define ROOTFIND_SOLNHIGHERTHANXMAX -4
|
||||
//@{
|
||||
//@}
|
||||
|
||||
|
||||
|
||||
|
||||
//! Root finder for 1D problems
|
||||
/*!
|
||||
*
|
||||
* The root finder solves a single nonlinear equation described below.
|
||||
*
|
||||
* \f[
|
||||
* f(x) = f_0
|
||||
* \f]
|
||||
*
|
||||
* \f$ f(x) \f$ is assumed to be single valued as a function of x.\f$ f(x) \f$ is not assumed to be continuous nor is
|
||||
* its derivative assumed to be well formed.
|
||||
*
|
||||
* Root finders are significantly different in the sense that do not have to rely
|
||||
* solely on Newton's method to find the answer to the problem. Instead they use a method to bound
|
||||
* the solution between high and low values and then use a method to refine that bound. The eventual
|
||||
* solution to the problem is presented as x_best and as a bound, delta_X, on the solution
|
||||
* component. Because of this, they are far more stable for functions and Jacobians that have discontinuities
|
||||
* or noise associated with them.
|
||||
*
|
||||
* The algorithm is a convolution of a local Secant method with an approach of finding a straddle in x.
|
||||
* The Jacobian is never required.
|
||||
*
|
||||
* There is a general breakdown of the algorithm into stages. The first stage seeks to find a straddle of the
|
||||
* function. The second stage seeks to reduce the bounds in x and f in order to satisfy the specification of the
|
||||
* stopping criteria. In the last stage the algorithm seeks to find the base value of x that satisfies the
|
||||
* original equation given what it current knows about the function.
|
||||
*
|
||||
* Globalization strategy
|
||||
*
|
||||
* Specifying the General Changes in x
|
||||
*
|
||||
* Supplying Hints with General Function Behavior Flags
|
||||
*
|
||||
*
|
||||
*
|
||||
* Stopping Criteria
|
||||
*
|
||||
* Specification of the Stopping Criteria
|
||||
*
|
||||
*
|
||||
* Additional constraints
|
||||
*
|
||||
* Bounds Criteria For the Routine
|
||||
*
|
||||
* Example
|
||||
*
|
||||
* @code
|
||||
* // Define a residual. The definition of a residual involves a lot more work than is shown here.
|
||||
* ResidEval * ec;
|
||||
* // Instantiate the root finder with the residual to be solved, ec.
|
||||
* RootFind rf(&ec);
|
||||
* // Set the relative and absolute tolerancess for f and x.
|
||||
* rf.setTol(1.0E-5, 1.0E-10, 1.0E-5, 1.0E-11);
|
||||
* // Give a hint about the function's dependence on x. This is needed, for example, if the function has
|
||||
* // flat regions.
|
||||
* rf.setFuncIsGenerallyIncreasing(true);
|
||||
* rf.setDeltaX(0.01);
|
||||
* // Supply an initial guess for the solution
|
||||
* double xbest = phiM;
|
||||
* double oldP = printLvl_;
|
||||
* // Set the print level for the solver. Zero produces no output. Two produces a summary table of each iteration.
|
||||
* rf.setPrintLvl(2);
|
||||
* // Define a minimum and maximum for the independent variable.
|
||||
* double phimin = 1.3;
|
||||
* double phimax = 2.2;
|
||||
* // Define a maximum iteration number
|
||||
* int itmax = 100;
|
||||
* // Define the f_0 value, and on return will contain the actual value of f(x) obtained
|
||||
* double currentObtained;
|
||||
* // Call the solver
|
||||
* status = rf.solve(phimin, phimax, 100, currentObtained, &xbest);
|
||||
* if (status == 0) {
|
||||
* if (printLvl_ > 1) {
|
||||
* printf("Electrode::integrateConstantCurrent(): Volts (%g amps) = %g\n", currentObtained, xbest);
|
||||
* }
|
||||
* } else {
|
||||
* if (printLvl_) {
|
||||
* printf("Electrode::integrateConstantCurrent(): bad status = %d Volts (%g amps) = %g\n",
|
||||
* status, currentObtained, xbest);
|
||||
* }
|
||||
* }
|
||||
* @endcode
|
||||
*
|
||||
* @todo Noise
|
||||
* @todo General Search to be done when all else fails
|
||||
*
|
||||
*/
|
||||
class RootFind {
|
||||
|
|
@ -110,14 +202,18 @@ namespace Cantera {
|
|||
//! 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
|
||||
* are close to one another. This is defined as being within factor * delXMeaningful() of each other.
|
||||
*
|
||||
* The basic premise here is that if the two numbers are too close, the noise
|
||||
* will prevent an accurate calculation of the function and its slope.
|
||||
*
|
||||
* @param x1 First number
|
||||
* @param x2 second number
|
||||
* @param factor Multiplicative factor to multiple deltaX with
|
||||
*
|
||||
* @return Returns a boolean indicating whether the two numbers are the same or not.
|
||||
*/
|
||||
bool theSame(doublereal x2, doublereal x1) const;
|
||||
bool theSame(doublereal x2, doublereal x1, doublereal factor = 1.0) const;
|
||||
|
||||
public:
|
||||
|
||||
|
|
@ -233,6 +329,9 @@ namespace Cantera {
|
|||
*/
|
||||
void setDeltaXMax(doublereal deltaX);
|
||||
|
||||
//! Print the iteration history table
|
||||
void printTable();
|
||||
|
||||
public:
|
||||
|
||||
//! Pointer to the residual function evaluator
|
||||
|
|
@ -249,6 +348,7 @@ namespace Cantera {
|
|||
|
||||
//! Relative tolerance for the value of f and x
|
||||
doublereal m_rtolf;
|
||||
|
||||
//! Relative tolerance for the value of x
|
||||
doublereal m_rtolx;
|
||||
|
||||
|
|
@ -256,10 +356,19 @@ namespace Cantera {
|
|||
doublereal m_maxstep;
|
||||
|
||||
protected:
|
||||
|
||||
//! Print level
|
||||
/*!
|
||||
* 0 No printing of any kind
|
||||
* 1 Single print line indicating success or failure of the routine.
|
||||
* 2 Summary table printed at the end of the routine, with a convergence history
|
||||
* 3 Printouts during the iteration are added. Summary table is printed out at the end.
|
||||
* if writeLogAllowed_ is turned on, a file is written out with the convergence history.
|
||||
*/
|
||||
int printLvl;
|
||||
|
||||
public:
|
||||
|
||||
//! Boolean to turn on the possibility of writing a log file.
|
||||
bool writeLogAllowed_;
|
||||
|
||||
|
|
@ -304,6 +413,57 @@ namespace Cantera {
|
|||
//! Internal variable tracking f(x) of smallest x tried.
|
||||
doublereal fx_minTried_;
|
||||
|
||||
|
||||
//! Structure containing the iteration history
|
||||
struct rfTable {
|
||||
//@{
|
||||
int its;
|
||||
int TP_its;
|
||||
double slope;
|
||||
double xval;
|
||||
double fval;
|
||||
int foundPos;
|
||||
int foundNeg;
|
||||
double deltaXConverged;
|
||||
double deltaFConverged;
|
||||
double delX;
|
||||
|
||||
|
||||
std::string reasoning;
|
||||
|
||||
void clear() {
|
||||
its = 0;
|
||||
TP_its = 0;
|
||||
slope = -1.0E300;
|
||||
xval = -1.0E300;
|
||||
fval = -1.0E300;
|
||||
reasoning = "";
|
||||
};
|
||||
|
||||
rfTable() :
|
||||
its(-2),
|
||||
TP_its(0),
|
||||
slope(-1.0E300),
|
||||
xval(-1.0E300),
|
||||
fval(-1.0E300),
|
||||
foundPos(0),
|
||||
foundNeg(0),
|
||||
deltaXConverged(-1.0E300),
|
||||
deltaFConverged(-1.0E300),
|
||||
delX(-1.0E300),
|
||||
reasoning("")
|
||||
{
|
||||
};
|
||||
|
||||
|
||||
//@}
|
||||
};
|
||||
|
||||
//! Vector of iteration histories
|
||||
std::vector<struct rfTable> rfHistory_;
|
||||
|
||||
|
||||
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue