From aa4dcf393c7c953950daed33260f753bd8df97c1 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Mon, 24 Mar 2014 04:03:44 +0000
Subject: [PATCH] [1D] Avoid Jacobian update after successful solve

Calling solve() again after a successful solution no longer automatically
triggers a Jacobian update. This enables relatively efficient sensitivity
analysis by sequentially perturbing the reaction multipliers and re-solving
the system. Since the perturbed system is close to the orignal, the solution
can be found after only a few steps, even when using the original Jacobian.
---
 include/cantera/oneD/StFlow.h | 22 ++++++++++++++++++++--
 src/oneD/MultiNewton.cpp      |  2 +-
 src/oneD/StFlow.cpp           |  4 ++--
 3 files changed, 23 insertions(+), 5 deletions(-)

diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h
index 8f58ebade..f838511ab 100644
--- a/include/cantera/oneD/StFlow.h
+++ b/include/cantera/oneD/StFlow.h
@@ -183,31 +183,49 @@ public:
     }
 
     void solveEnergyEqn(size_t j=npos) {
+        bool changed = false;
         if (j == npos)
             for (size_t i = 0; i < m_points; i++) {
+                if (!m_do_energy[i]) {
+                    changed = true;
+                }
                 m_do_energy[i] = true;
             }
         else {
+            if (!m_do_energy[j]) {
+                changed = true;
+            }
             m_do_energy[j] = true;
         }
         m_refiner->setActive(0, true);
         m_refiner->setActive(1, true);
         m_refiner->setActive(2, true);
-        needJacUpdate();
+        if (changed) {
+            needJacUpdate();
+        }
     }
 
     void fixTemperature(size_t j=npos) {
+        bool changed = false;
         if (j == npos)
             for (size_t i = 0; i < m_points; i++) {
+                if (m_do_energy[i]) {
+                    changed = true;
+                }
                 m_do_energy[i] = false;
             }
         else {
+            if (m_do_energy[j]) {
+                changed = true;
+            }
             m_do_energy[j] = false;
         }
         m_refiner->setActive(0, false);
         m_refiner->setActive(1, false);
         m_refiner->setActive(2, false);
-        needJacUpdate();
+        if (changed) {
+            needJacUpdate();
+        }
     }
 
     bool doSpecies(size_t k) {
diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp
index d06e222a9..a9b07006b 100644
--- a/src/oneD/MultiNewton.cpp
+++ b/src/oneD/MultiNewton.cpp
@@ -382,7 +382,6 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
     int nJacReeval = 0;
 
     while (1 > 0) {
-
         // Check whether the Jacobian should be re-evaluated.
         if (jac.age() > m_maxAge) {
             writelog("\nMaximum Jacobian age reached ("+int2str(m_maxAge)+")\n", loglevel);
@@ -427,6 +426,7 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
 
         // convergence
         else if (m == 1) {
+            jac.setAge(0); // for efficient sensitivity analysis
             break;
         }
 
diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
index 6d9d3cb30..20e8d471d 100644
--- a/src/oneD/StFlow.cpp
+++ b/src/oneD/StFlow.cpp
@@ -209,11 +209,11 @@ void StFlow::_finalize(const doublereal* x)
     bool e = m_do_energy[0];
     for (j = 0; j < m_points; j++) {
         if (e || nz == 0) {
-            setTemperature(j, T(x, j));
+            m_fixedtemp[j] = T(x, j);
         } else {
             zz = (z(j) - z(0))/(z(m_points - 1) - z(0));
             tt = linearInterp(zz, m_zfix, m_tfix);
-            setTemperature(j, tt);
+            m_fixedtemp[j] = tt;
         }
         for (k = 0; k < m_nsp; k++) {
             setMassFraction(j, k, Y(x, k, j));