diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 2c9ff4cd9..9cc67f1ee 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -108,7 +108,6 @@ namespace Cantera { m_numTotalLinearSolves(0), m_numTotalNewtIts(0), m_min_newt_its(0), - filterNewstep(0), maxNewtIts_(100), m_jacFormMethod(NSOLN_JAC_NUM), m_nJacEval(0), @@ -178,7 +177,6 @@ namespace Cantera { m_numTotalLinearSolves(0), m_numTotalNewtIts(0), m_min_newt_its(0), - filterNewstep(0), maxNewtIts_(100), m_jacFormMethod(NSOLN_JAC_NUM), m_nJacEval(0), @@ -237,7 +235,6 @@ namespace Cantera { m_numTotalLinearSolves = right.m_numTotalLinearSolves; 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; @@ -1261,6 +1258,12 @@ namespace Cantera { m_normSolnFRaw = solnErrorNorm(DATA_PTR(stp), "Initial Undamped Step of the iteration", 0); } calcSolnToResNormVector(); + + /* + * Filter out bad directions + */ + doublereal normFilter = filterNewStep(time_curr, DATA_PTR(m_y_n), DATA_PTR(stp)); + // Damp the Newton step /* @@ -1328,7 +1331,7 @@ namespace Cantera { bool m_filterIntermediate = false; if (m_filterIntermediate) { if (m == 0) { - (void) filterNewStep(time_n, DATA_PTR(y_new), DATA_PTR(ydot_new)); + (void) filterNewSolution(time_n, DATA_PTR(y_new), DATA_PTR(ydot_new)); } } @@ -1592,7 +1595,7 @@ namespace Cantera { ydot[j] += dy * CJ; } /* - * Call the functon + * Call the function */ @@ -1618,17 +1621,17 @@ namespace Cantera { if (m_print_flag >= 7 && s_print_NumJac) { if (neq_ < 30) { - printf("\t\tCurrent Matrix:\n"); + printf("\t\tCurrent Matrix and Residual:\n"); printf("\t\t I,J | "); for (j = 0; j < neq_; j++) { printf(" %5d ", j); } - printf("\n"); + printf("| Residual \n"); printf("\t\t --"); for (j = 0; j < neq_; j++) { printf("------------"); } - printf("\n"); + printf("| -----------\n"); for (i = 0; i < neq_; i++) { @@ -1636,14 +1639,14 @@ namespace Cantera { for (j = 0; j < neq_; j++) { printf(" % 11.4E", J(i,j) ); } - printf("\n"); + printf(" | % 11.4E\n", f[i]); } printf("\t\t --"); for (j = 0; j < neq_; j++) { printf("------------"); } - printf("\n"); + printf("--------------\n"); } } @@ -1681,7 +1684,7 @@ namespace Cantera { } } //==================================================================================================================== - // Apply a filtering step + // Apply a filtering process to the new step /* * @param timeCurrent Current value of the time * @param y_current current value of the solution @@ -1689,8 +1692,24 @@ namespace Cantera { * * @return Returns the norm of the value of the amount filtered */ - doublereal NonlinearSolver::filterNewStep(const doublereal timeCurrent, doublereal * const y_current, doublereal *const ydot_current) { - return 0.0; + doublereal NonlinearSolver::filterNewStep(const doublereal timeCurrent, + const doublereal * const ybase, doublereal * const step0) { + doublereal tmp = m_func->filterNewStep(timeCurrent, ybase, step0); + return tmp; + } + //==================================================================================================================== + // Apply a filtering process to the new solution + /* + * @param timeCurrent Current value of the time + * @param y_current current value of the solution + * @param ydot_current Current value of the solution derivative. + * + * @return Returns the norm of the value of the amount filtered + */ + doublereal NonlinearSolver::filterNewSolution(const doublereal timeCurrent, + doublereal * const y_current, doublereal *const ydot_current) { + doublereal tmp = m_func->filterSolnPrediction(timeCurrent, y_current); + return tmp; } //==================================================================================================================== // Compute the Residual Weights diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index a09245728..65136e4c7 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -247,7 +247,16 @@ namespace Cantera { doublereal * const ydot, int num_newt_its); - //! Apply a filtering step + //! Apply a filtering process to the step + /*! + * @param timeCurrent Current value of the time + * @param ybase current value of the solution + * + * @return Returns the norm of the value of the amount filtered + */ + doublereal filterNewStep(const doublereal timeCurrent, const doublereal * const ybase, doublereal * const step0); + + //! Apply a filter to the solution /*! * @param timeCurrent Current value of the time * @param y_current current value of the solution @@ -255,8 +264,8 @@ namespace Cantera { * * @return Returns the norm of the value of the amount filtered */ - doublereal filterNewStep(const doublereal timeCurrent, doublereal * const y_current, doublereal * const ydot_current); - + doublereal filterNewSolution(const doublereal timeCurrent, doublereal * const y_current, doublereal * const ydot_current); + //! Return the factor by which the undamped Newton step 'step0' //! must be multiplied in order to keep the update within the bounds of an accurate jacobian. @@ -543,9 +552,6 @@ namespace Cantera { //! Minimum number of newton iterations to use int m_min_newt_its; - //! Boolean that turns on solution filtering - int filterNewstep; - //! Maximum number of newton iterations int maxNewtIts_; diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp index 220a64560..e27aa8db9 100644 --- a/Cantera/src/numerics/ResidJacEval.cpp +++ b/Cantera/src/numerics/ResidJacEval.cpp @@ -210,6 +210,20 @@ namespace Cantera { } //==================================================================================================================== // Filter the solution predictions + /* + * Codes might provide a predicted step change. This routine filters the predicted + * solution vector eliminating illegal directions. + * + * @param t Time (input) + * @param y Solution vector (input, output) + * @param step Proposed step in the solution that will be cropped + */ + doublereal ResidJacEval::filterNewStep(doublereal t, const doublereal * const ybase, doublereal * const step) + { + return 0.0; + } + //==================================================================================================================== + // Filter the solution predictions /* * Codes might provide a predicted solution vector. This routine filters the predicted * solution vector. @@ -217,8 +231,9 @@ namespace Cantera { * @param t Time (input) * @param y Solution vector (input, output) */ - void ResidJacEval::filterSolnPrediction(doublereal t, doublereal * const y) + doublereal ResidJacEval::filterSolnPrediction(doublereal t, doublereal * const y) { + return 0.0; } //==================================================================================================================== // Evalulate any stopping criteria other than a final time limit diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h index ee4c2a7d0..7ddb4e0a9 100644 --- a/Cantera/src/numerics/ResidJacEval.h +++ b/Cantera/src/numerics/ResidJacEval.h @@ -132,6 +132,18 @@ namespace Cantera { */ virtual void getInitialConditions(const doublereal t0, doublereal * const y, doublereal * const ydot); + //! Filter the solution predictions + /*! + * Codes might provide a predicted step change. This routine filters the predicted + * solution vector eliminating illegal directions. + * + * @param t Time (input) + * @param y Solution vector (input, output) + * @param step Proposed step in the solution that will be cropped + */ + virtual doublereal filterNewStep(const doublereal t, const doublereal * const ybase, + doublereal * const step); + //! Filter the solution predictions /*! * Codes might provide a predicted solution vector. This routine filters the predicted @@ -140,7 +152,7 @@ namespace Cantera { * @param t Time (input) * @param y Solution vector (input, output) */ - virtual void filterSolnPrediction(const doublereal t, doublereal * const y); + virtual doublereal filterSolnPrediction(const doublereal t, doublereal * const y); //! Set a global value of the absolute tolerance /*!