From f727cd44af203f74348afcc179f835a963804690 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 5 Jan 2011 03:55:13 +0000 Subject: [PATCH] Added in a more robust treatment of handling limits --- Cantera/src/numerics/NonlinearSolver.cpp | 20 +++++++++++--- Cantera/src/numerics/RootFind.cpp | 35 ++++++++++++++++++++++-- Cantera/src/numerics/RootFind.h | 6 ++-- 3 files changed, 52 insertions(+), 9 deletions(-) diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 1d7df013a..602892d68 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -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]); diff --git a/Cantera/src/numerics/RootFind.cpp b/Cantera/src/numerics/RootFind.cpp index b72ffefff..e6e703ede 100644 --- a/Cantera/src/numerics/RootFind.cpp +++ b/Cantera/src/numerics/RootFind.cpp @@ -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; } //================================================================================================ diff --git a/Cantera/src/numerics/RootFind.h b/Cantera/src/numerics/RootFind.h index 98e3a3ede..5e6dd6429 100644 --- a/Cantera/src/numerics/RootFind.h +++ b/Cantera/src/numerics/RootFind.h @@ -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