From 6e8f90c59542c16f987be0d3623a867cda5aa5b5 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Tue, 16 May 2017 12:52:31 +0900
Subject: [PATCH] working but little slow

---
 EEqn.H          |  2 ++
 YEqn.H          |  2 ++
 createFields.H  | 40 +++++++++++++++++++++++++++++++
 eReactingFoam.C |  2 --
 neEqn.H         | 63 ++++++++++++++++++++++++++++++-------------------
 5 files changed, 83 insertions(+), 26 deletions(-)

diff --git a/EEqn.H b/EEqn.H
index 9262b91..431f4ac 100644
--- a/EEqn.H
+++ b/EEqn.H
@@ -33,4 +33,6 @@
 
     Info<< "min/max(T) = "
         << min(T).value() << ", " << max(T).value() << endl;
+
+    updateElectronTransport = true;
 }
diff --git a/YEqn.H b/YEqn.H
index ac1e073..fe2bb39 100644
--- a/YEqn.H
+++ b/YEqn.H
@@ -178,4 +178,6 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
+
+    updateElectronTransport = true;
 }
diff --git a/createFields.H b/createFields.H
index 1bfe560..85d39a1 100644
--- a/createFields.H
+++ b/createFields.H
@@ -1,3 +1,5 @@
+Switch updateElectronTransport(true);
+
 Info<< "Reading physicalProperties\n" << endl;
 
 IOdictionary physicalProperties
@@ -546,3 +548,41 @@ forAll (neutrals, iidx)
     neutralFluxBFs.set(iidx,
             new surfaceScalarField::GeometricBoundaryField (phi.boundaryField()));
 }
+
+
+Info<< "Reading physicalProperties/crossSections\n" << endl;
+dictionary crossSections
+(
+    physicalProperties.subDict("crossSections")
+);
+
+wordList csSpecies((crossSections.lookup("species")));
+Info<< "Cross sections are considered for following species\n" << endl;
+Info<< csSpecies << endl;
+
+labelList csSpeciesI(csSpecies.size(), -1);
+
+PtrList< DataEntry< scalar > > csList (csSpecies.size());
+
+forAll (csSpecies, isp)
+{
+    csSpeciesI[isp] = composition.species()[csSpecies[isp]];
+
+    csList.set(
+        isp,
+        DataEntry<scalar>::New(csSpecies[isp], crossSections)
+    );
+}
+
+Info<< csSpeciesI << endl;
+
+scalarField csSpeciesW (csSpecies.size(), 0.0);
+forAll (csSpecies, isp)
+{
+    csSpeciesW[isp] = composition.W(csSpeciesI[isp]);
+}
+
+Info<< csSpeciesW << endl;
+
+
+scalarField sigma_m (Te);
diff --git a/eReactingFoam.C b/eReactingFoam.C
index e5a6cdd..0952458 100644
--- a/eReactingFoam.C
+++ b/eReactingFoam.C
@@ -39,8 +39,6 @@ Description
 
 #include "CSV.H"
 
-// #include "bolos.h"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
diff --git a/neEqn.H b/neEqn.H
index dd1978f..6dcc746 100644
--- a/neEqn.H
+++ b/neEqn.H
@@ -1,35 +1,49 @@
 {
     // Electron swarm parameter
 
-    EnTd = En.internalField();
-    EnTd *= EnToTableUnit;
-
-    Te.internalField() = TeOfEn.value(EnTd) * TeFac;
-    forAll(rho, celli)
+    if (updateElectronTransport)
     {
-        Te[celli] = max(Te[celli], T[celli]);
-    }
-    Te.correctBoundaryConditions();
+        Te = T;
+        Te.correctBoundaryConditions();
 
-    if (mobility_f_of_Te)
-    {
-        EnTd = Te.internalField();
-        EnTd *= TeToTableUnit;
-    }
+        scalarField X (csSpecies.size(), 0.0);
+        scalarField &W (csSpeciesW);
 
-    mue.internalField() = mueN.value(EnTd) * mueNFac;
+        forAll (Te, cid)
+        {
+            scalar eEps = (3./2.) * kB.value() * Te[cid]
+                / eCharge.value();
+
+            sigma_m[cid] = 0;
+            forAll (csSpecies, isp)
+            {
+                X[isp] = Y[csSpeciesI[isp]][cid]/W[isp];
+                sigma_m[cid] +=
+                    X[isp] * csList[isp].value(eEps);
+            }
+            sigma_m[cid] /= sum(X);
+        }
+
+        Info <<  min(sigma_m) << " / " << max(sigma_m) << endl;
+
+/*
+        mue.internalField() = eCharge.value()
+            / (sqrt(eMass.value() * kB.value() * Te) * sigma_m);
+        mue.correctBoundaryConditions();
+*/
+        mue.internalField() = mueNFac * eCharge.value() / 3.0
+            / (sqrt(pi * eMass.value() * kB.value() * Te / 8.0) * sigma_m);
+        mue.correctBoundaryConditions();
 
-    if (calculateDe)
-    {
         De = mue * Te * (kB / eCharge);
-    }
-    else
-    {
-        De.internalField() = DeN.value(EnTd) * DeNFac;
-    }
 
-    mue.correctBoundaryConditions();
-    De.correctBoundaryConditions();
+        De.correctBoundaryConditions();
+
+        updateElectronTransport = false;
+
+        mue.writeMinMax(Info);
+        De.writeMinMax(Info);
+    }
 
 
     Udrift = -(mue/ng)*E;
@@ -37,7 +51,8 @@
     {
         Udrift *= relaxDrift;
     }
-    Uthermal = -((De/ng/Te)*fvc::grad(Te));
+    // Uthermal = -((De/ng/Te)*fvc::grad(Te));
+    Uthermal = -(fvc::grad(De/ng));
     ve = (linearInterpolate(Udrift+Uthermal) & mesh.Sf()) + q;
 
     const surfaceScalarField &msf = mesh.magSf();