From 09dba5d5d556fcaa4ddf033317c189243d179f49 Mon Sep 17 00:00:00 2001
From: tlichtenegger <thomas.lichtenegger@jku.at>
Date: Mon, 9 Apr 2018 12:30:31 +0200
Subject: [PATCH] Additional options for treatment of Y fields.

---
 .../solvers/cfdemSolverRhoPimpleChem/YEqn.H   | 29 +++++++++++--------
 .../cfdemSolverRhoPimpleChem/createFields.H   |  6 +++-
 .../cfdemSolverRhoPimpleChem/debugYEqn.H      |  2 +-
 3 files changed, 23 insertions(+), 14 deletions(-)

diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H b/applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H
index b8af296d..3c82ddc8 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H
@@ -19,10 +19,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     forAll(Y, i)
     {
-        if (Y[i].name() != inertSpecie || ignoreInertSpecie)
+        if (Y[i].name() == inertSpecie) inertIndex = i;
+        if (Y[i].name() != inertSpecie || propagateInertSpecie)
         {
             volScalarField& Yi = Y[i];
-	    
+
             fvScalarMatrix YiEqn
             (
                 fvm::ddt(rhoeps, Yi)
@@ -42,23 +43,27 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
             fvOptions.correct(Yi);
 
-             #include "debugYEqn.H"
-
             Yi.max(0.0);
-            Yt += Yi;
-        }
-        else
-        {
-            inertIndex = i;
+            if (Y[i].name() != inertSpecie) Yt += Yi;
         }
     }
 
-    if (ignoreInertSpecie)
+    if (inertIndex!=-1)
     {
+        Y[inertIndex].max(inertLowerBound);
+        Y[inertIndex].min(inertUpperBound);
+    }
+
+    if (propagateInertSpecie)
+    {
+        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + VSMALL);
         forAll(Y,i)
         {
-            volScalarField& Yi = Y[i];
-            Yi = Yi/Yt;
+            if (i!=inertIndex)
+            {
+                volScalarField& Yi = Y[i];
+                Yi = Yi/(Yt+VSMALL);
+            }
         }
     }
     else
diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/createFields.H b/applications/solvers/cfdemSolverRhoPimpleChem/createFields.H
index edf31d39..477e00d6 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/createFields.H
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/createFields.H
@@ -17,10 +17,14 @@
     // read molecular weight
     volScalarField W(composition.W());
 
-    bool ignoreInertSpecie = true;
+    bool propagateInertSpecie = true;
 
     const word inertSpecie(thermo.lookup("inertSpecie"));
 
+    const scalar inertLowerBound(thermo.lookupOrDefault<scalar>("inertLowerBound",0.0));
+
+    const scalar inertUpperBound(thermo.lookupOrDefault<scalar>("inertUpperBound",1.0));
+
     if (!composition.contains(inertSpecie))
     {
         FatalErrorIn(args.executable())
diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/debugYEqn.H b/applications/solvers/cfdemSolverRhoPimpleChem/debugYEqn.H
index 2a475fd7..4e6dc253 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/debugYEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/debugYEqn.H
@@ -21,7 +21,7 @@
             Info << "\nartificial mass of species " << Y[i].name() << " per time step: "<< fvc::domainIntegrate(artMass) << endl;
             if(lVCell > -1)
             {
-                Pout << "lowest value of " << lowestValue << " at cell " << lVCell << " with coordinates" << endl;
+                Pout << Y[i].name() << ": time / lowest value " << runTime.timeName() << "\t" << lowestValue << "\n\tat cell " << lVCell << " with coordinates";
                 Pout << "\t" << mesh.C()[lVCell].component(0) << "\t" << mesh.C()[lVCell].component(1)  << "\t" << mesh.C()[lVCell].component(2) << endl;
             }
 }