From 0c178dbdd8576d65ced462c917a4e36d47d1e404 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 10 Feb 2011 20:59:53 +0000 Subject: [PATCH] Prototyping a Dogleg method. Added in the calcualtion of the steepest descent and the Cauchy point --- Cantera/src/numerics/NonlinearSolver.cpp | 98 +++++++++++++++++++++++- Cantera/src/numerics/NonlinearSolver.h | 21 +++++ 2 files changed, 118 insertions(+), 1 deletion(-) diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index b632f72d1..b57c7cb36 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -122,6 +122,11 @@ namespace Cantera { atolk_(0), m_print_flag(0), m_ScaleSolnNormToResNorm(0.001) +#ifdef DEBUG_DOGLEG + ,descentDir_(0), + residNorm2Cauchy_(0.0), + Jd_(0) +#endif { neq_ = m_func->nEquations(); @@ -145,6 +150,13 @@ namespace Cantera { atolk_[i] = atolBase_; m_ewt[i] = atolk_[i]; } + +#ifdef DEBUG_DOGLEG + jacCopy_.resize(neq_, neq_, 0.0); + descentDir_.resize(neq_, 0.0); + Jd_.resize(neq_, 0.0); +#endif + } //==================================================================================================================== NonlinearSolver::NonlinearSolver(const NonlinearSolver &right) : @@ -191,6 +203,11 @@ namespace Cantera { atolk_(0), m_print_flag(0), m_ScaleSolnNormToResNorm(0.001) +#ifdef DEBUG_DOGLEG + ,descentDir_(0), + residNorm2Cauchy_(0.0), + Jd_(0) +#endif { *this =operator=(right); } @@ -248,6 +265,11 @@ namespace Cantera { atolk_ = right.atolk_; m_print_flag = right.m_print_flag; m_ScaleSolnNormToResNorm = right.m_ScaleSolnNormToResNorm; +#ifdef DEBUG_DOGLEG + jacCopy_ = right.jacCopy_; + descentDir_ = right.descentDir_; + Jd_ = right.Jd_; +#endif return *this; } @@ -704,6 +726,75 @@ namespace Cantera { return info; } //==================================================================================================================== + #ifdef DEBUG_DOGLEG + // Do a steepest descent calculation + /* + * This call must be made on the unfactored jacobian! + */ + int NonlinearSolver::doCauchyPointSolve(SquareMatrix& jac) + { + double rowFac = 1.0; + // Calculate desDir = -0.5 * R dot J + /* + * this would be faster:: + * vector_fp &dd = jac.data(); + * descentDir[j] -= 0.5 * resid[i] * dd[i*neq_ * j[; + */ + for (int j = 0; j < neq_; j++) { + descentDir_[j] = 0.0; + double colFac = 1.0; + if (m_colScaling) { + colFac = 1.0/m_colScales[j]; + } + for (int i = 0; i < neq_; i++) { + if (m_rowScaling) { + rowFac = 1.0/m_rowScales[i]; + } + descentDir_[j] -= 0.5 * m_resid[i] * jac.value(i,j) *colFac / (m_residWts[i] * m_residWts[i]); + } + } + for (int j = 0; j < neq_; j++) { + Jd_[j] = 0.0; + double colFac = 1.0; + if (m_colScaling) { + colFac = 1.0/m_colScales[j]; + } + for (int i = 0; i < neq_; i++) { + if (m_rowScaling) { + rowFac = 1.0/m_rowScales[i]; + } + Jd_[j] += descentDir_[j] * jac.value(i,j) *rowFac * colFac/ m_residWts[i]; + } + } + double RJd_norm = 0.0; + double JdJd_norm = 0.0; + for (int i = 0; i < neq_; i++) { + RJd_norm += m_resid[i] * Jd_[i] / m_residWts[i]; + JdJd_norm += Jd_[i] * Jd_[i]; + } + double lambda = - RJd_norm / (JdJd_norm); + + for (int i = 0; i < neq_; i++) { + descentDir_[i] *= lambda; + } + + residNorm2Cauchy_ = m_normResidFRaw * m_normResidFRaw - RJd_norm * RJd_norm / (JdJd_norm*JdJd_norm); + + // Compute the weighted norm of the undamped step size descentDir_[] + doublereal sDD = solnErrorNorm(DATA_PTR(descentDir_), "SteepestDescentDir", 10); + + + if (m_print_flag > 2) { + printf("\t\t\tdoCauchyPointSolve: Steepest descent to Cauchy point: \n"); + printf("\t\t\t Rraw = %g Rpred = %g, deltaX = %g\n", m_normResidFRaw, residNorm2Cauchy_, sDD); + } + + + + return 0; + } +#endif + //==================================================================================================================== void NonlinearSolver::setDefaultDeltaBoundsMagnitudes() { for (int i = 0; i < neq_; i++) { @@ -1174,7 +1265,9 @@ namespace Cantera { int m = 0; bool forceNewJac = false; doublereal s1=1.e30; - +#ifdef DEBUG_DOGLEG + jacCopy_ = jac; +#endif // std::vector y_curr(neq_, 0.0); std::vector ydot_curr(neq_, 0.0); @@ -1298,6 +1391,9 @@ namespace Cantera { m_normResid0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 0, DATA_PTR(m_y_n)); } +#ifdef DEBUG_DOGLEG + doCauchyPointSolve(jac); +#endif // compute the undamped Newton step info = doNewtonSolve(time_curr, DATA_PTR(m_y_n), DATA_PTR(ydot_curr), DATA_PTR(stp), jac, m_print_flag); diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index aeec870e5..5f93f5aab 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -496,6 +496,11 @@ namespace Cantera { //! solution norms. void calcSolnToResNormVector(); +#ifdef DEBUG_DOGLEG + int doCauchyPointSolve(SquareMatrix& jac); + +#endif + //! Set the print level from the rootfinder /*! * @@ -702,6 +707,22 @@ namespace Cantera { //! Scale factor for turning residual norms into solution norms double m_ScaleSolnNormToResNorm; + +#ifdef DEBUG_DOGLEG + //! Copy of the jacobian that doesn't get overwritten when the inverse is determined + SquareMatrix jacCopy_; + + //! Steepest descent direction. This is also the distance to the Cauchy Point + std::vector descentDir_; + + //! Expected value of the residual norm at the Cauchy point + doublereal residNorm2Cauchy_; + + //! Jacobian times the Steepest descent direction. + std::vector Jd_; + +#endif + public: //! Turn off printing of time /*!