From bc492823fcf8d047548a59f7a64c8a3bc11cb9bc Mon Sep 17 00:00:00 2001
From: tlichtenegger <thomas.lichtenegger@jku.at>
Date: Tue, 20 Oct 2020 12:01:25 +0200
Subject: [PATCH] Restructuring the solver.

---
 .../rcfdemSolverRhoSteadyPimpleChem/UEqn.H    | 40 +++++++-------
 .../rcfdemSolverRhoSteadyPimpleChem/YEqn.H    |  6 +--
 .../createFields.H                            |  2 +
 .../rcfdemSolverRhoSteadyPimpleChem/pEqn.H    |  5 ++
 .../rcfdemSolverRhoSteadyPimpleChem.C         | 52 +++++++++----------
 5 files changed, 57 insertions(+), 48 deletions(-)

diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H
index 9e75193e..44cb1f6e 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H
@@ -1,7 +1,7 @@
 // Solve the Momentum equation
 particleCloud.otherForces(fOther);
 
-tmp<fvVectorMatrix> tUEqn
+fvVectorMatrix UEqn
 (
     fvm::div(phi, U)
   + particleCloud.divVoidfractionTau(U, voidfractionRec)
@@ -10,24 +10,26 @@ tmp<fvVectorMatrix> tUEqn
  ==
     fvOptions(rho, U)
 );
-fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn.relax();
-
-fvOptions.constrain(UEqn);
-
-if (modelType=="B" || modelType=="Bfull")
+if (stepcounter%nEveryFlow==0)
 {
-    solve(UEqn == -fvc::grad(p)+ Ksl*UsRec);
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (modelType=="B" || modelType=="Bfull")
+    {
+        solve(UEqn == -fvc::grad(p)+ Ksl*UsRec);
+    }
+    else
+    {
+        solve(UEqn == -voidfractionRec*fvc::grad(p)+ Ksl*UsRec);
+    }
+
+
+    #include "limitU.H"
+
+    fvOptions.correct(U);
+
+    K = 0.5*magSqr(U);
 }
-else
-{
-    solve(UEqn == -voidfractionRec*fvc::grad(p)+ Ksl*UsRec);
-}
-
-//U.relax();
-#include "limitU.H"
-
-fvOptions.correct(U);
-
-K = 0.5*magSqr(U);
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H
index de5c56e9..c5cbaca0 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H
@@ -30,13 +30,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
             fvScalarMatrix YiEqn
             (
-                //fvm::ddt(rhoeps, Yi)
-              //+
                 mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(voidfractionRec*turbulence->muEff(), Yi)
               ==
                 combustion->R(Yi)
-              + particleCloud.chemistryM(0).Smi(i)
+              + particleCloud.chemistryM(0).Smi(i)*p/p.prevIter()
               + fvOptions(rho, Yi)
             );
 
@@ -46,6 +44,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
             YiEqn.solve(mesh.solver("Yi"));
 
+            Yi.relax();
+
             fvOptions.correct(Yi);
 
             Yi.max(0.0);
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H
index c0ac42dc..1556801e 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H
@@ -212,6 +212,8 @@ Info<< "Reading thermophysical properties\n" << endl;
          linearInterpolate(rho*U*voidfraction) & mesh.Sf()
     );
 
+    bool transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
+
     dimensionedScalar rhoMax
     (
         dimensionedScalar::lookupOrDefault
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H
index 6b40829c..c8e4c697 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H
@@ -3,6 +3,9 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
+if (stepcounter%nEveryFlow==0)
+{
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
 if (modelType=="A")
@@ -89,3 +92,5 @@ else
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
+
+}
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C
index 9e22b7f8..e2dc841c 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C
@@ -95,7 +95,7 @@ int main(int argc, char *argv[])
 
     const IOdictionary& couplingProps = particleCloud.couplingProperties();
     label nEveryFlow(couplingProps.lookupOrDefault<label>("nEveryFlow",1));
-    Info << "Solving flow equations every " << nEveryFlow << " steps.\n" << endl;
+    Info << "Solving flow equations for U and p every " << nEveryFlow << " steps.\n" << endl;
     label stepcounter = 0;
 
     Info<< "\nStarting time loop\n" << endl;
@@ -103,6 +103,7 @@ int main(int argc, char *argv[])
     scalar m(0.0);
     scalar m0(0.0);
     label counter(0);
+    p.storePrevIter();
 
     while (runTime.run())
     {
@@ -144,39 +145,38 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().start(26,"Flow");
         volScalarField rhoeps("rhoeps",rho*voidfractionRec);
-        if (stepcounter%nEveryFlow==0)
+
+        while (pimple.loop())
         {
-            while (pimple.loop())
-            {
-                // if needed, perform drag update here
+            // if needed, perform drag update here
 #if OPENFOAM_VERSION_MAJOR < 6
-                if (pimple.nCorrPIMPLE() <= 1)
+            if (pimple.nCorrPIMPLE() <= 1)
 #else
-                if (pimple.nCorrPimple() <= 1)
+            if (pimple.nCorrPimple() <= 1)
 #endif
-                {
-                    #include "rhoEqn.H"
-                }
+            {
+                #include "rhoEqn.H"
+            }
 
-                // --- Pressure-velocity PIMPLE corrector loop
+            // --- Pressure-velocity PIMPLE corrector loop
 
-                #include "UEqn.H"
-                #include "YEqn.H"
-                #include "EEqn.H"
 
-                // --- Pressure corrector loop
-                while (pimple.correct())
-                {
-                    // besides this pEqn, OF offers a "pimple consistent"-option
-                    #include "molConc.H"
-                    #include "pEqn.H"
-                    rhoeps=rho*voidfractionRec;
-                }
+            #include "UEqn.H"
+            #include "EEqn.H"
 
-                if (pimple.turbCorr())
-                {
-                    turbulence->correct();
-                }
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                // besides this pEqn, OF offers a "pimple consistent"-option
+                #include "molConc.H"
+                #include "pEqn.H"
+                rhoeps=rho*voidfractionRec;
+            }
+            #include "YEqn.H"
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
             }
         }