From 86b85b3eefd0ea6d700127fef7889cff44a86b94 Mon Sep 17 00:00:00 2001
From: Yeongdo Park <ignis@postech.ac.kr>
Date: Tue, 1 Jan 2019 14:54:54 -0500
Subject: [PATCH] binary diffusivity and pure viscosity test

---
 .../diffusivityModel/diffusivityModel.C       | 56 +++++++++++++++++++
 .../diffusivityModel/diffusivityModel.H       |  4 ++
 testApp/chemFoam.C                            | 39 ++++++++++++-
 3 files changed, 97 insertions(+), 2 deletions(-)

diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index d840384..83f0dd9 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -714,6 +714,62 @@ void Foam::diffusivityModel::correct()
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::diffusivityModel::Dbinary
+(const scalar T, const scalar p, const UList<scalar> &Y)
+{
+    const speciesTable &species_(thermo_.composition().species());
+
+    scalarSymmetricSquareMatrix Dij(species_.size());
+
+    scalarField localY(Y);
+
+    scalarField muI(species_.size());
+
+    {
+        const scalar pi = p;
+        const scalar Ti = T;
+
+        GasState state (pi, Ti, localY);
+
+        calculateMuD (muI, Dij, state);
+
+    }
+
+    tmp<scalarField> pDii(new scalarField(UList<scalar>(Dij.v(), Dij.size())));
+
+    return pDii;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::diffusivityModel::mupure(const scalar T, const scalar p, const UList<scalar> &Y)
+{
+    const speciesTable &species_(thermo_.composition().species());
+
+    scalarSymmetricSquareMatrix Dij(species_.size());
+
+    scalarField localY(Y);
+
+    scalarField Di(species_.size());
+    scalarField Dii(species_.size());
+    scalarField *pmuI = new scalarField(species_.size());
+    scalarField &muI = *pmuI;
+
+    {
+        const scalar pi = p;
+        const scalar Ti = T;
+
+        GasState state (pi, Ti, localY);
+
+        const scalarField &localX = state.X();
+
+        calculateMuD (muI, Dij, state);
+
+    }
+
+    return tmp<scalarField>(pmuI);
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::diffusivityModel::kpure(const scalar T, const scalar p, const UList<scalar> &Y)
 {
     const speciesTable &species_(thermo_.composition().species());
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index 0c4ca8d..7eec9e8 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.H
@@ -217,6 +217,10 @@ public:
 
         // Check
 
+        tmp<scalarField> Dbinary(const scalar T, const scalar p, const UList<scalar> &Y);
+
+        tmp<scalarField> mupure(const scalar T, const scalar p, const UList<scalar> &Y);
+
         tmp<scalarField> kpure(const scalar T, const scalar p, const UList<scalar> &Y);
 
         // Edit
diff --git a/testApp/chemFoam.C b/testApp/chemFoam.C
index 048bf9b..8c7a8b8 100644
--- a/testApp/chemFoam.C
+++ b/testApp/chemFoam.C
@@ -176,6 +176,8 @@ int main(int argc, char *argv[])
         thermo.correct();
         diff.correct();
 
+        scalarField DijO (diff.Dbinary(thermo.T()[0], p[0], XY));
+
         // Cantera gas transport
         gas_.setState_TPY(thermo.T()[0], p[0], XY.data());
 
@@ -196,14 +198,14 @@ int main(int argc, char *argv[])
             }
         }
 
+        tr_->getBinaryDiffCoeffs(nCanteraSp_, Dij.data());
+
         if (idxFracOne < 0)
         {
             tr_->getMixDiffCoeffsMass(Dc.data());
         }
         else
         {
-            tr_->getBinaryDiffCoeffs(nCanteraSp_, Dij.data());
-
             forAll(Dc, k)
             {
                 Dc[k] = Dij[nCanteraSp_*idxFracOne+k];
@@ -213,6 +215,11 @@ int main(int argc, char *argv[])
 
         scalarField kSp (diff.kpure(thermo.T()[0], p[0], XY));
 
+        scalarField muSpO (diff.mupure(thermo.T()[0], p[0], XY));
+
+        scalarField muSpC (kSp);
+
+        tr_->getSpeciesViscosities(muSpC.data());
 
 
         forAll(thermo.composition().species(), k)
@@ -227,6 +234,22 @@ int main(int argc, char *argv[])
                 ;
         }
 
+        forAll(thermo.composition().species(), k)
+        {
+        forAll(thermo.composition().species(), j)
+        {
+            logWriter
+                (
+                 thermo.composition().species()[k] + "=>" +
+                 thermo.composition().species()[j],
+                 thermo.T()[0],
+                 DijO[nCanteraSp_*k+j],
+                 Dij[nCanteraSp_*k+j]
+                )
+                ;
+        }
+        }
+
         logWriter
             (
              "mu",
@@ -236,6 +259,18 @@ int main(int argc, char *argv[])
             )
             ;
 
+        forAll(thermo.composition().species(), k)
+        {
+            logWriter
+                (
+                 "mu_"+thermo.composition().species()[k],
+                 thermo.T()[0],
+                 muSpO[k],
+                 muSpC[k]
+                )
+                ;
+        }
+
         logWriter
             (
              "k",