diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 596d44cda..1d7df013a 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -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; + } + //===================================================================================================================== } diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index 08699226e..96dd0e08a 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -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 & lowBoundsConstraintVector(); - //! return an editable vector of the high bounds constraints + //! Return an editable vector of the high bounds constraints std::vector & 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 diff --git a/Cantera/src/numerics/ResidEval.h b/Cantera/src/numerics/ResidEval.h index 1c24843e7..da65c4cfa 100644 --- a/Cantera/src/numerics/ResidEval.h +++ b/Cantera/src/numerics/ResidEval.h @@ -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 diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp index b1ae4edf1..cb6d08988 100644 --- a/Cantera/src/numerics/ResidJacEval.cpp +++ b/Cantera/src/numerics/ResidJacEval.cpp @@ -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; } //==================================================================================================================== diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h index e1b1fc069..08a23fc35 100644 --- a/Cantera/src/numerics/ResidJacEval.h +++ b/Cantera/src/numerics/ResidJacEval.h @@ -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: diff --git a/Cantera/src/numerics/RootFind.cpp b/Cantera/src/numerics/RootFind.cpp index 0bc16bd9d..b72ffefff 100644 --- a/Cantera/src/numerics/RootFind.cpp +++ b/Cantera/src/numerics/RootFind.cpp @@ -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); diff --git a/Cantera/src/numerics/RootFind.h b/Cantera/src/numerics/RootFind.h index c8168ed89..98e3a3ede 100644 --- a/Cantera/src/numerics/RootFind.h +++ b/Cantera/src/numerics/RootFind.h @@ -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_; }; }