From 4d245501cf8dffb35a92882f22f62ffde681fd43 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Mon, 30 May 2011 17:57:53 +0000
Subject: [PATCH] Update on the NonlinearSolver:     Update to
 :calcSolnToResNormVector() to align it with documentation.

---
 Cantera/src/numerics/NonlinearSolver.cpp | 51 +++++++++++++++++-------
 Cantera/src/numerics/NonlinearSolver.h   | 14 +++++--
 2 files changed, 46 insertions(+), 19 deletions(-)

diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp
index 9190e9911..467d82138 100644
--- a/Cantera/src/numerics/NonlinearSolver.cpp
+++ b/Cantera/src/numerics/NonlinearSolver.cpp
@@ -348,17 +348,12 @@ namespace Cantera {
    *  The program always assumes that atol is specific
    *  to the solution component
    *
-   * param y  vector of the current solution values
+   * @param y  vector of the current solution values
    */
   void NonlinearSolver::createSolnWeights(const doublereal * const y) {
     for (int i = 0; i < neq_; i++) {
       m_ewt[i] = rtol_ * fabs(y[i]) + atolk_[i];
     }
-#ifdef DEBUG_DOGLEG
-    //  for (int i = 0; i < neq_; i++) {
-    //   m_ewt[i] = 1.0E-4;
-    // }
-#endif
   }
   //====================================================================================================================
   // set bounds constraints for all variables in the problem
@@ -665,8 +660,10 @@ namespace Cantera {
       }
       for (jcol = 0; jcol < neq_; jcol++) {
 	for (irow = 0; irow < neq_; irow++) {
-	  m_rowScales[irow] += fabs(*jptr);
+	  //m_rowScales[irow] += fabs(*jptr);
 	  if (m_colScaling) {
+	    // This is needed in order to mitgate the change in J_ij carried out just above this loop.
+	    // Alternatively, we could move this loop up to the top
 	    m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol] / m_colScales[jcol];
 	  } else {
 	    m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol];
@@ -699,16 +696,40 @@ namespace Cantera {
 
   }
   //====================================================================================================================
+  // Calculate the scaling factor for translating residual norms into solution norms.
+  /*
+   *  This routine calls computeResidWts() a couple of times in the calculation of m_ScaleSolnNormToResNorm.
+   *  A more sophisticated routine may do more with signs to get a better value. Perhaps, a series of calculations
+   *  with different signs attached may be in order. Then, m_ScaleSolnNormToResNorm would be calculated
+   *  as the minimum of a series of calculations.
+   */
   void NonlinearSolver::calcSolnToResNormVector()
   { 
-    //  double oldVal =  m_ScaleSolnNormToResNorm;
-    // if (m_normSolnFRaw > 1.0E-13) {
-    // m_ScaleSolnNormToResNorm = 1.0E-2 * m_normResid0 / m_normSolnFRaw * oldVal;
-    //} 
-    // m_normResid0 = m_normSolnFRaw;
-    //computeResidWts();
-    m_ScaleSolnNormToResNorm = 1.0;
-    //double tmp = residErrorNorm(DATA_PTR(m_resid));
+    if (! jacCopy_.m_factored) { 
+      m_ScaleSolnNormToResNorm = 1.0;
+      computeResidWts();
+      for (int n = 0; n < neq_; n++) {
+        m_wksp[n] = 0.0;
+      }   
+      doublereal *jptr = &(*(jacCopy_.begin()));
+      for (int jcol = 0; jcol < neq_; jcol++) {
+	for (int irow = 0; irow < neq_; irow++) {
+	  m_wksp[irow] += (*jptr) * m_ewt[jcol];
+	  jptr++;
+	}
+      }
+      double resNormOld = residErrorNorm(DATA_PTR(m_wksp));
+      if (resNormOld > 0.0) {
+	m_ScaleSolnNormToResNorm = m_ScaleSolnNormToResNorm * resNormOld;
+      }
+      if (m_ScaleSolnNormToResNorm < 1.0E-8) {
+	m_ScaleSolnNormToResNorm = 1.0E-8;
+      }
+      // Recalculate the residual weights now that we know the value of m_ScaleSolnNormToResNorm 
+      computeResidWts();
+    } else {
+      throw CanteraError("NonlinearSolver::calcSolnToResNormVector()" , "Logic error");
+    }
   }
   //====================================================================================================================
   // Compute the undamped Newton step based on the current jacobian and an input rhs
diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h
index 5ea6ddb6d..454286edf 100644
--- a/Cantera/src/numerics/NonlinearSolver.h
+++ b/Cantera/src/numerics/NonlinearSolver.h
@@ -579,8 +579,13 @@ namespace Cantera {
      */
     void setMaxNewtIts(const int maxNewtIts);
 
-    //! Calculate the scaling factor for translating residual norms into 
-    //! solution norms.
+    //! Calculate the scaling factor for translating residual norms into solution norms.
+    /*!
+     *  This routine calls computeResidWts() a couple of times in the calculation of m_ScaleSolnNormToResNorm.
+     *  A more sophisticated routine may do more with signs to get a better value. Perhaps, a series of calculations
+     *  with different signs attached may be in order. Then, m_ScaleSolnNormToResNorm would be calculated
+     *  as the minimum of a series of calculations.
+     */
     void calcSolnToResNormVector();
 
     //! Calculate the steepest descent direction and the Cauchy Point where the quadratic formulation 
@@ -758,8 +763,9 @@ namespace Cantera {
 
     //! Weights for normalizing the values of the residuals
     /*!
-     *  These are computed if row scaling, m_rowScaling, is turned on. They are calculated currently as the
-     *  sum of the absolute values  jacobian multiplied by the solution weight function
+     *  They are calculated as the  sum of the absolute values of the jacobian 
+     *  multiplied by the solution weight function.
+     *  This is carried out in scaleMatrix().
      */
     std::vector<doublereal> m_rowWtScales;