From 48526118f22112066e688bcdd1d38f2e992b787e Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Wed, 18 Feb 2009 22:31:09 +0000
Subject: [PATCH] Added more functionality to the stefan-maxwel solve.

---
 Cantera/src/transport/LiquidTransport.cpp | 58 +++++++++++++++++++++++
 Cantera/src/transport/LiquidTransport.h   |  6 +++
 2 files changed, 64 insertions(+)

diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp
index 45afd1f62..6882415b8 100644
--- a/Cantera/src/transport/LiquidTransport.cpp
+++ b/Cantera/src/transport/LiquidTransport.cpp
@@ -336,6 +336,7 @@ namespace Cantera {
    for (int i = 0; i < itop; i++) {
      m_Grad_X[i] = grad_X[i];
    }
+   update_Grad_lnAC();
  }
 
 
@@ -595,6 +596,63 @@ namespace Cantera {
   }
 
 
+  // We formulate the directional derivative
+  void LiquidTransport::update_Grad_lnAC() {
+    int k;
+    
+
+    for (int a = 0; a < m_nDim; a++) {
+      // We form the directional derivative
+      double * ma_Grad_X = &m_Grad_X[a*m_nsp];
+      double sum = 0.0;
+      for (k = 0; k < m_nsp; k++) {
+        sum += ma_Grad_X[k] * ma_Grad_X[k];
+      }
+      if (sum == 0.0) {
+	for (k = 0; k < m_nsp; k++) {
+	  m_Grad_lnAC[m_nsp * a + k] = 0.0;
+	}
+	continue;
+      }
+      double mag = 1.0E-7 / sum;
+    
+	for (k = 0; k < m_nsp; k++) {
+	  Xdelta_[k] = m_molefracs[k] + mag * ma_Grad_X[k];
+	  if (Xdelta_[k] > 1.0) {
+	    Xdelta_[k] = 1.0;
+	  }
+	  if (Xdelta_[k] < 0.0) {
+	    Xdelta_[k] = 0.0;
+	  }
+	}
+      m_thermo->setMoleFractions(DATA_PTR(Xdelta_));
+      m_thermo->getActivityCoefficients(DATA_PTR(lnActCoeffMolarDelta_));
+      for (k = 0; k < m_nsp; k++) {
+	lnActCoeffMolarDelta_[k] = log(lnActCoeffMolarDelta_[k]);
+      }
+
+      for (k = 0; k < m_nsp; k++) {
+	m_Grad_lnAC[m_nsp * a + k] = sum * (lnActCoeffMolarDelta_[k] - log(actCoeffMolar_[k])) / mag;
+      }
+    }
+    m_thermo->setMoleFractions(DATA_PTR(m_molefracs));
+    double Tbase = m_thermo->temperature();
+    double T_new = Tbase - 1.0E-6;
+    m_thermo->setTemperature(T_new);
+    
+     m_thermo->getActivityCoefficients(DATA_PTR(lnActCoeffMolarDelta_));
+     double *dlnActCoeffdT =  &Xdelta_[0];
+     for (k = 0; k < m_nsp; k++) {
+       dlnActCoeffdT[k] ==  (lnActCoeffMolarDelta_[k] - log(actCoeffMolar_[k]))/(-1.0E-6);
+     }
+     for (int a = 0; a < m_nDim; a++) {
+       for (k = 0; k < m_nsp; k++) {
+	 m_Grad_lnAC[m_nsp * a + k] +=  dlnActCoeffdT[k] * m_Grad_T[a];
+       }
+     }
+    m_thermo->setTemperature(Tbase);
+  }
+
   /*************************************************************************
    *
    *    methods to update temperature-dependent properties
diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h
index 372ff7cdd..2d7b4f4eb 100644
--- a/Cantera/src/transport/LiquidTransport.h
+++ b/Cantera/src/transport/LiquidTransport.h
@@ -254,6 +254,8 @@ namespace Cantera {
      */
     virtual void set_Grad_X(const doublereal* const grad_X);
 
+   virtual void update_Grad_lnAC();
+
   protected:
     //! Handles the effects of changes in the Temperature, internally
     //! within the object.
@@ -499,6 +501,8 @@ namespace Cantera {
      */
     vector_fp m_molefracs;
 
+    vector_fp Xdelta_;
+
     //! Local copy of the concentrations of the species in the phase
     /*!
      * Update info?
@@ -521,6 +525,8 @@ namespace Cantera {
 
     vector_fp actCoeffMolar_;
 
+    vector_fp lnActCoeffMolarDelta_;
+
     //! Stefan-Maxwell Diffusion Coefficients at T, P and C
     /*!
      *   These diffusion coefficients are considered to be