From 92f8fc42928e1d7d409cb4f73f493beed51d28ef Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Sun, 7 Sep 2014 21:17:03 +0000
Subject: [PATCH] [1D] Change multicomponent diffusion calculation to improve
 convergence

Change the formulation of multicomponent diffusion fluxes to avoid using the
equivalent Fickian diffusion coefficient. That formulation leads to negative
diffusivities for some species at some grid points, which generate numerical
instabilities while timestepping. This formulation produces a more dense
Jacobian, but with a much better eigenvalue spectrum.
---
 src/oneD/StFlow.cpp | 48 ++++++++++++++++++++++-----------------------
 1 file changed, 24 insertions(+), 24 deletions(-)

diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
index 41d975f6a..56c621ebb 100644
--- a/src/oneD/StFlow.cpp
+++ b/src/oneD/StFlow.cpp
@@ -415,34 +415,21 @@ void StFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
             m_tcon[j] = m_trans->thermalConductivity();
         }
     } else if (m_transport_option == c_Multi_Transport) {
-        doublereal sum, sumx, wtm, dz;
-        doublereal eps = 1.0e-12;
-        for (size_t m = j0; m < j1; m++) {
-            setGasAtMidpoint(x,m);
-            dz = m_z[m+1] - m_z[m];
-            wtm = m_thermo->meanMolecularWeight();
-
-            m_visc[m] = (m_dovisc ? m_trans->viscosity() : 0.0);
-
-            m_trans->getMultiDiffCoeffs(m_nsp,
-                                        DATA_PTR(m_multidiff) + mindex(0,0,m));
+        for (size_t j = j0; j < j1; j++) {
+            setGasAtMidpoint(x,j);
+            doublereal wtm = m_thermo->meanMolecularWeight();
+            doublereal rho = m_thermo->density();
+            m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
+            m_trans->getMultiDiffCoeffs(m_nsp, &m_multidiff[mindex(0,0,j)]);
 
+            // Use m_diff as storage for the factor outside the summation
             for (size_t k = 0; k < m_nsp; k++) {
-                sum = 0.0;
-                sumx = 0.0;
-                for (size_t j = 0; j < m_nsp; j++) {
-                    if (j != k) {
-                        sum += m_wt[j]*m_multidiff[mindex(k,j,m)]*
-                               ((X(x,j,m+1) - X(x,j,m))/dz + eps);
-                        sumx += (X(x,j,m+1) - X(x,j,m))/dz;
-                    }
-                }
-                m_diff[k + m*m_nsp] = sum/(wtm*(sumx+eps));
+                m_diff[k+j*m_nsp] = m_wt[k] * rho / (wtm*wtm);
             }
 
-            m_tcon[m] = m_trans->thermalConductivity();
+            m_tcon[j] = m_trans->thermalConductivity();
             if (m_do_soret) {
-                m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + m*m_nsp);
+                m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + j*m_nsp);
             }
         }
     }
@@ -506,7 +493,6 @@ void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
     switch (m_transport_option) {
 
     case c_Mixav_Transport:
-    case c_Multi_Transport:
         for (j = j0; j < j1; j++) {
             sum = 0.0;
             wtm = m_wtm[j];
@@ -525,6 +511,20 @@ void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
         }
         break;
 
+    case c_Multi_Transport:
+        for (j = j0; j < j1; j++) {
+            dz = z(j+1) - z(j);
+
+            for (k = 0; k < m_nsp; k++) {
+                doublereal sum = 0.0;
+                for (size_t m = 0; m < m_nsp; m++) {
+                    sum += m_wt[m] * m_multidiff[mindex(k,m,j)] * (X(x,m,j+1)-X(x,m,j));
+                }
+                m_flux(k,j) = sum * m_diff[k+j*m_nsp] / dz;
+            }
+        }
+        break;
+
     default:
         throw CanteraError("updateDiffFluxes","unknown transport model");
     }