Updated to OF 4.x.

applications/solvers/cfdemSolverIB/Make/options

@@ -1,24 +1,27 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude
+
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lfvOptions \

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -39,7 +39,8 @@ Contributions
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
 
 #include "cfdemCloudIB.H"
 #include "implicitCouple.H"
@@ -52,10 +53,9 @@ Contributions
 
 #include "cellSet.H"
 
-#if defined(version22)
-    #include "meshToMeshNew.H"
-    #include "fvIOoptionList.H"
-#endif
+#include "meshToMeshNew.H"
+#include "fvIOoptionList.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,14 +66,14 @@ int main(int argc, char *argv[])
     #include "createTime.H"
 
     #include "createDynamicFvMesh.H"
+  
+    #include "createControl.H"
 
     #include "createFields.H"
 
     #include "initContinuityErrs.H"
 
-    #if defined(version22)
     #include "createFvOptions.H"
-    #endif
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
@@ -91,7 +91,6 @@ int main(int argc, char *argv[])
         interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
         mesh.update(); //dyM
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
@@ -107,45 +106,33 @@ int main(int argc, char *argv[])
                 fvm::ddt(voidfraction,U)
               + fvm::div(phi, U)
               + turbulence->divDevReff(U)
-                #if defined(version22)
                 ==
                 fvOptions(U)
-                #endif
             );
 
             UEqn.relax();
 
-            #if defined(version22)
             fvOptions.constrain(UEqn);
-            #endif
 
-            if (momentumPredictor)
+            if (piso.momentumPredictor())
             {
                 solve(UEqn == -fvc::grad(p));
             }
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+           while (piso.correct())
             {
                 volScalarField rUA = 1.0/UEqn.A();
                 surfaceScalarField rUAf(fvc::interpolate(rUA));
 
                 U = rUA*UEqn.H();
-                #ifdef version23
+
                 phi = (fvc::interpolate(U) & mesh.Sf())
                     + rUAf*fvc::ddtCorr(U, phi);
-                #else
-                phi = (fvc::interpolate(U) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rUA, U, phi);
-                #endif
+
                 adjustPhi(phi, U, p);
 
-                #if defined(version22)
-                fvOptions.relativeFlux(phi);
-                #endif
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correct())
                 {
                     // Pressure corrector
 
@@ -156,20 +143,9 @@ int main(int argc, char *argv[])
 
                     pEqn.setReference(pRefCell, pRefValue);
 
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
+                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (piso.finalNonOrthogonalIter())
                     {
                         phi -= pEqn.flux();
                     }

applications/solvers/cfdemSolverPiso/Make/options

@@ -1,19 +1,22 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C

@@ -36,7 +36,9 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
@@ -51,7 +53,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     // create cfdemCloud
@@ -67,7 +71,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
@@ -85,7 +88,7 @@ int main(int argc, char *argv[])
 
        //Force Checks
        vector fTotal(0,0,0);
-       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField()));
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
        reduce(fImpTotal, sumOp<vector>());
        Info << "TotalForceExp: " << fTotal << endl;
        Info << "TotalForceImp: " << fImpTotal << endl;
@@ -100,95 +103,19 @@ int main(int argc, char *argv[])
             // Pressure-velocity PISO corrector
             {
                 // Momentum predictor
-                fvVectorMatrix UEqn
-                (
-                    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
-                  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
-//                + turbulence->divDevReff(U)
-                  + particleCloud.divVoidfractionTau(U, voidfraction)
-                  ==
-                  - fvm::Sp(Ksl/rho,U)
-                );
-
-                UEqn.relax();
-                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
-                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else if (momentumPredictor)
-                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+                 #include "UEqn.H"
 
                 // --- PISO loop
 
-                //for (int corr=0; corr<nCorr; corr++)
-                int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-                for (int corr=0; corr<nCorrSoph; corr++)
+                while (piso.correct())
                 {
-                    volScalarField rUA = 1.0/UEqn.A();
-
-                    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                    volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
-                    surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
-
-                    U = rUA*UEqn.H();
-
-                    #ifdef version23
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + rUAfvoidfraction*fvc::ddtCorr(U, phi);
-                    #else
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                    #endif
-                    surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                    surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
-                    if (modelType=="A")
-                        rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                    // Non-orthogonal pressure corrector loop
-                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                    {
-                        // Pressure corrector
-                        fvScalarMatrix pEqn
-                        (
-                            fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                        );
-                        pEqn.setReference(pRefCell, pRefValue);
-
-                        if
-                        (
-                            corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                        {
-                            pEqn.solve(mesh.solver("pFinal"));
-                        }
-                        else
-                        {
-                            pEqn.solve();
-                        }
-
-                        if (nonOrth == nNonOrthCorr)
-                        {
-                            phiGes -= pEqn.flux();
-                            phi = phiGes;
-                        }
-
-                    } // end non-orthogonal corrector loop
-
-                    #include "continuityErrorPhiPU.H"
-
-                    if (modelType=="B" || modelType=="Bfull")
-                        U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                    else
-                        U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                    U.correctBoundaryConditions();
-
-                } // end piso loop
+                    #include "pEqn.H"
+                }
             }
 
+            laminarTransport.correct();
             turbulence->correct();
-        }// end solveFlow
+        }
         else
         {
             Info << "skipping flow solution." << endl;

applications/solvers/cfdemSolverPiso/createFields.H

@@ -122,3 +122,5 @@ surfaceScalarField phi
     (
         incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
+
+#include "createMRF.H"

applications/solvers/cfdemSolverPisoMS/Make/options

@@ -1,16 +1,23 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I ../cfdemSolverPiso \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -l$(CFDEM_LIB_NAME)
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C

@@ -25,23 +25,26 @@ License
     along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cfdemSolverPisoMS
+    cfdemSolverPiso
 
 Description
     Transient solver for incompressible flow.
     Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
-    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, 
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
     where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloudMS.H"
 #include "implicitCouple.H"
 #include "clockModel.H"
 #include "smoothingModel.H"
+#include "forceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +53,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     // create cfdemCloud
@@ -59,118 +64,62 @@ int main(int argc, char *argv[])
     #include "checkModelType.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+    Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
-        Info<< "\nStarting time loop\n" << endl;
-            particleCloud.clockM().start(1,"Global");
+        particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
         particleCloud.clockM().start(2,"Coupling");
-        particleCloud.evolve(voidfraction,Us,U);
-      
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+    
         Info << "update Ksl.internalField()" << endl;
         Ksl = particleCloud.momCoupleM(0).impMomSource();
-        particleCloud.smoothingM().smoothen(Ksl);
         Ksl.correctBoundaryConditions();
 
+       //Force Checks
+       vector fTotal(0,0,0);
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+       reduce(fImpTotal, sumOp<vector>());
+       Info << "TotalForceExp: " << fTotal << endl;
+       Info << "TotalForceImp: " << fImpTotal << endl;
+
         #include "solverDebugInfo.H"
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
-        // Pressure-velocity PISO corrector
+
+        if(particleCloud.solveFlow())
         {
-            // Momentum predictor
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(voidfraction,U) //particleCloud.ddtVoidfractionU(U,voidfraction) //
-              + fvm::div(phi, U)
-//              + turbulence->divDevReff(U)
-              + particleCloud.divVoidfractionTau(U, voidfraction)
-             ==
-              - fvm::Sp(Ksl/rho,U)
-            );
-
-            if (modelType=="B")
-                UEqn == - fvc::grad(p) + Ksl/rho*Us;
-            else
-                UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
-
-            UEqn.relax();
-
-            if (momentumPredictor)
-                solve(UEqn);
-
-            // --- PISO loop
-
-            //for (int corr=0; corr<nCorr; corr++)
-            int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-            for (int corr=0; corr<nCorrSoph; corr++)
+            // Pressure-velocity PISO corrector
             {
-                volScalarField rUA = 1.0/UEqn.A();
+                // Momentum predictor
+                 #include "UEqn.H"
 
-                surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+                // --- PISO loop
 
-                U = rUA*UEqn.H();
-
-                phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() );
-                     //+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
-                if (modelType=="A")
-                    rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correct())
                 {
-                    // Pressure corrector
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                    );
-                    pEqn.setReference(pRefCell, pRefValue);
+                    #include "pEqn.H"
+                }
+            }
 
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
-
-                    if (nonOrth == nNonOrthCorr)
-                    {
-                        phiGes -= pEqn.flux();
-                    }
-
-                } // end non-orthogonal corrector loop
-
-                #include "continuityErrorPhiPU.H"
-
-                if (modelType=="B")
-                    U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                else
-                    U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                U.correctBoundaryConditions();
-
-            } // end piso loop
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
         }
-
-        turbulence->correct();
 
         runTime.write();
 
@@ -183,7 +132,7 @@ int main(int argc, char *argv[])
     }
 
     Info<< "End\n" << endl;
-    
+
     return 0;
 }
 

applications/solvers/cfdemSolverPisoScalar/Make/options

@@ -1,19 +1,23 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I ../cfdemSolverPiso \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C

@@ -25,7 +25,7 @@ License
     along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cfdemSolverPisoScalar
+    cfdemSolverPiso
 
 Description
     Transient solver for incompressible flow.
@@ -36,10 +36,13 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
+#include "clockModel.H"
 #include "smoothingModel.H"
 #include "forceModel.H"
 
@@ -50,7 +53,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     // create cfdemCloud
@@ -62,12 +67,14 @@ int main(int argc, char *argv[])
     Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
+        particleCloud.clockM().start(1,"Global");
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
         if(hasEvolved)
@@ -79,119 +86,38 @@ int main(int argc, char *argv[])
         Ksl = particleCloud.momCoupleM(0).impMomSource();
         Ksl.correctBoundaryConditions();
 
+       //Force Checks
+       vector fTotal(0,0,0);
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+       reduce(fImpTotal, sumOp<vector>());
+       Info << "TotalForceExp: " << fTotal << endl;
+       Info << "TotalForceImp: " << fImpTotal << endl;
 
         #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
 
-        // get scalar source from DEM        
-        particleCloud.forceM(1).manipulateScalarField(Tsource);
-        Tsource.correctBoundaryConditions();
-
-        // solve scalar transport equation
-        fvScalarMatrix TEqn
-        (
-           fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
-         + fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
-         - fvm::laplacian(DT*voidfraction, T)
-         ==
-           Tsource
-        );
-        TEqn.relax();
-        TEqn.solve();
+        particleCloud.clockM().start(26,"Flow");
+	
+        #include "TEqn.H"
 
         if(particleCloud.solveFlow())
         {
             // Pressure-velocity PISO corrector
             {
                 // Momentum predictor
-                fvVectorMatrix UEqn
-                (
-                    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
-                  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
-//                + turbulence->divDevReff(U)
-                  + particleCloud.divVoidfractionTau(U, voidfraction)
-                  ==
-                  - fvm::Sp(Ksl/rho,U)
-                );
-
-                UEqn.relax();
-                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
-                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else if (momentumPredictor)
-                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+                 #include "UEqn.H"
 
                 // --- PISO loop
 
-                //for (int corr=0; corr<nCorr; corr++)
-                int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-                for (int corr=0; corr<nCorrSoph; corr++)
+                while (piso.correct())
                 {
-                    volScalarField rUA = 1.0/UEqn.A();
-
-                    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                    volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
-                    surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
-
-                    U = rUA*UEqn.H();
-
-                    #ifdef version23
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + rUAfvoidfraction*fvc::ddtCorr(U, phi);
-                    #else
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                    #endif
-                    surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                    surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
-                    if (modelType=="A")
-                        rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                    // Non-orthogonal pressure corrector loop
-                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                    {
-                        // Pressure corrector
-                        fvScalarMatrix pEqn
-                        (
-                            fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                        );
-                        pEqn.setReference(pRefCell, pRefValue);
-
-                        if
-                        (
-                            corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                        {
-                            pEqn.solve(mesh.solver("pFinal"));
-                        }
-                        else
-                        {
-                            pEqn.solve();
-                        }
-
-                        if (nonOrth == nNonOrthCorr)
-                        {
-                            phiGes -= pEqn.flux();
-                            phi = phiGes;
-                        }
-
-                    } // end non-orthogonal corrector loop
-
-                    #include "continuityErrorPhiPU.H"
-
-                    if (modelType=="B" || modelType=="Bfull")
-                        U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                    else
-                        U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                    U.correctBoundaryConditions();
-
-                } // end piso loop
+                    #include "pEqn.H"
+                }
             }
 
+            laminarTransport.correct();
             turbulence->correct();
-        }// end solveFlow
+        }
         else
         {
             Info << "skipping flow solution." << endl;
@@ -202,6 +128,9 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
     }
 
     Info<< "End\n" << endl;

applications/solvers/cfdemSolverPisoScalar/createFields.H

@@ -1,36 +1,36 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading physical velocity field U" << endl;
-    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-//========================
-// drag law modelling
-//========================
-
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading physical velocity field U" << endl;
+    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// drag law modelling
+//========================
+
     Info<< "\nReading momentum exchange field Ksl\n" << endl;
     volScalarField Ksl
     (
@@ -44,8 +44,8 @@
         ),
         mesh
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.0)
-    );
-
+    );
+
     Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (
@@ -58,8 +58,8 @@
             IOobject::AUTO_WRITE
         ),
         mesh
-    );
-
+    );
+
     Info<< "\nCreating density field rho\n" << endl;
     volScalarField rho
     (
@@ -71,27 +71,27 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh,
+        mesh,
         dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
-    );
-
-    Info<< "Reading particle velocity field Us\n" << endl;
-    volVectorField Us
-    (
-        IOobject
-        (
-            "Us",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-//========================
-// scalar field modelling
-//========================
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// scalar field modelling
+//========================
     Info<< "\nCreating dummy density field rho = 1\n" << endl;
     volScalarField T
     (
@@ -103,10 +103,10 @@
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh//,
+        mesh//,
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 273.15)
-    );
-
+    );
+
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Tsource
     (
@@ -118,57 +118,59 @@
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh//,
+        mesh//,
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 0.0)
-    );
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar DT
-    (
-        transportProperties.lookup("DT")
-    );
-
-//========================
-
-//#   include "createPhi.H"
-#ifndef createPhi_H
-#define createPhi_H
-Info<< "Reading/calculating face flux field phi\n" << endl;
-surfaceScalarField phi
- (
-     IOobject
-     (
-        "phi",
-        runTime.timeName(),
-         mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-     ),
-     linearInterpolate(U*voidfraction) & mesh.Sf()
- );
-#endif
-
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::turbulenceModel> turbulence
-    (
-        incompressible::turbulenceModel::New(U, phi, laminarTransport)
-    );
+    );
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar DT
+    (
+        transportProperties.lookup("DT")
+    );
+
+//========================
+
+//#   include "createPhi.H"
+#ifndef createPhi_H
+#define createPhi_H
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+ (
+     IOobject
+     (
+        "phi",
+        runTime.timeName(),
+         mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+     ),
+     linearInterpolate(U*voidfraction) & mesh.Sf()
+ );
+#endif
+
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+#include "createMRF.H"

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -51,4 +51,4 @@
     fvOptions.correct(he);
 
     thermo.correct();
-}
+}

applications/solvers/cfdemSolverRhoPimple/Make/options

@@ -4,23 +4,24 @@ PFLAGS+= -Dcompre
 
 EXE_INC = \
      $(PFLAGS) \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
+    -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \
-    -lcompressibleTurbulenceModel \
-    -lcompressibleRASModels \
-    -lcompressibleLESModels \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
     -lfiniteVolume \
     -lmeshTools \
     -lsampling \

applications/solvers/cfdemSolverRhoPimple/UEqn.H

@@ -1,34 +1,33 @@
 // Solve the Momentum equation
 particleCloud.otherForces(fOther);
 
-tmp<fvVectorMatrix> UEqn
-
+tmp<fvVectorMatrix> tUEqn
 (
-    fvm::ddt(rhoeps,U)
-   + fvm::div(phi, U)
-   + particleCloud.divVoidfractionTau(U, voidfraction)
-   + fvm::Sp(Ksl,U)
-   - fOther
+    fvm::ddt(rhoeps, U)
+  + fvm::div(phi, U)
+  + particleCloud.divVoidfractionTau(U, voidfraction)
+  + fvm::Sp(Ksl,U)
+  - fOther
  ==
     fvOptions(rho, U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
 
-fvOptions.constrain(UEqn());
 if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
 {
-    solve(UEqn() == -fvc::grad(p)+ Ksl*Us);
+    solve(UEqn == -fvc::grad(p)+ Ksl*Us);
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
 }
 else if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -voidfraction*fvc::grad(p)+ Ksl*Us);
+    solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us);
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
 }
-
-

applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C

@@ -23,18 +23,19 @@ Description
     Transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
     algorithm.
     Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
-    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 2.3,
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
     where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
-
 #include "fvCFD.H"
 #include "psiThermo.H"
-#include "turbulenceModel.H"
+#include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
-#include "fvIOoptionList.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
 
 #include "cfdemCloudEnergy.H"
 #include "implicitCouple.H"
@@ -49,24 +50,26 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "postProcess.H"
+  
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-  
-    bool UfromPhi(runTime.controlDict().lookupOrDefault<bool>("UfromPhi",false));
-    Info << "UfromPhi = " << UfromPhi << endl;
-
-    pimpleControl pimple(mesh);
-
-    #include "createFields.H"
-    #include "createPcorrTypes.H"
-    #include "createFvOptions.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
-    
+    #include "createFields.H"
+    #include "createFvOptions.H"
+
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
     cfdemCloudEnergy particleCloud(mesh);
     #include "checkModelType.H"
+    
+    turbulence->validate();
+  //        #include "compressibleCourantNo.H"
+  //  #include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -97,12 +100,12 @@ int main(int argc, char *argv[])
         Ksl = particleCloud.momCoupleM(0).impMomSource();
         Ksl.correctBoundaryConditions();
 
-       //Force Checks
-       vector fTotal(0,0,0);
-       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField()));
-       reduce(fImpTotal, sumOp<vector>());
-       Info << "TotalForceExp: " << fTotal << endl;
-       Info << "TotalForceImp: " << fImpTotal << endl;
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef()));
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
 
         #include "solverDebugInfo.H"
         particleCloud.clockM().stop("Coupling");
@@ -113,8 +116,8 @@ int main(int argc, char *argv[])
         {
             #include "rhoEqn.H"
         }
-        //volScalarField rhoeps("rhoeps",rho*voidfraction);
-	rhoeps=rho*voidfraction;
+        
+        volScalarField rhoeps("rhoeps",rho*voidfraction);
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
@@ -124,8 +127,9 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
+	        // besides this pEqn, OF offers a "pimple consistent"-option
                 #include "pEqn.H"
-	      	rhoeps=rho*voidfraction;
+	      	rhoeps=rho*voidfraction;	
             }
 
             if (pimple.turbCorr())

applications/solvers/cfdemSolverRhoPimple/createFields.H

@@ -114,8 +114,6 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
     );
 
-    #ifndef compressibleCreatePhi_H
-    #define compressibleCreatePhi_H
     Info<< "Reading/calculating face flux field phi\n" << endl;
     surfaceScalarField phi
     (
@@ -123,16 +121,34 @@ Info<< "Reading thermophysical properties\n" << endl;
         (
             "phi",
             runTime.timeName(),
-             mesh,
+            mesh,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
          ),
          linearInterpolate(rho*U*voidfraction) & mesh.Sf()
     );
-    #endif
 
-    dimensionedScalar rhoMax(pimple.dict().lookup("rhoMax"));
-    dimensionedScalar rhoMin(pimple.dict().lookup("rhoMin"));
+    dimensionedScalar rhoMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMax",
+            pimple.dict(),
+            dimDensity,
+            GREAT
+        )
+    );
+
+    dimensionedScalar rhoMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMin",
+            pimple.dict(),
+            dimDensity,
+            0
+        )
+    );
 
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
@@ -145,6 +161,8 @@ Info<< "Reading thermophysical properties\n" << endl;
             thermo
         )
     );
+    
+    mesh.setFluxRequired(p.name());
 
     Info<< "Creating field dpdt\n" << endl;
     volScalarField dpdt
@@ -192,74 +210,4 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh
     );
 
-//===============================
-
-volScalarField divphi
-    (
-        IOobject
-        (
-            "divphi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 0.0)
-    );
-    
-    volScalarField divphi2
-    (
-        IOobject
-        (
-            "divphi2",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 0.0)
-    );
-
-volVectorField gP
-    (
-        IOobject
-        (
-            "gP",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedVector("0", dimensionSet(1, -2, -2, 0, 0), vector::zero)
-);
-    
-volVectorField gP2
-    (
-        IOobject
-        (
-            "gP2",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedVector("0", dimensionSet(1, -2, -2, 0, 0), vector::zero)
-);
-
-    volScalarField rhoeps
-    (
-        IOobject
-        (
-            "rhoeps",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 0.0)
-    );
+//===============================

applications/solvers/cfdemSolverRhoPimple/pEqn.H

@@ -3,57 +3,24 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*voidfraction*rAU));
-
-volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn().H();
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
+if (modelType=="A")
+{
+    rhorAUf *= fvc::interpolate(voidfraction);
+}
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
-
 if (pimple.transonic())
 {
 //  transonic version not implemented yet
-  
-/*    surfaceScalarField phid
-    (
-        "phid",
-        fvc::interpolate(psi)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-        )
-    );
-
-    fvOptions.makeRelative(fvc::interpolate(psi), phid);
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAUf, p)
-          ==
-            fvOptions(psi, p, rho.name())
-        );
-
-        fvOptions.constrain(pEqn);
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi == pEqn.flux();
-        }
-    }
-*/
 }
 else
 {
@@ -61,24 +28,16 @@ else
     (
         "phiHbyA",
         (
-            (fvc::interpolate(rhoeps*HbyA) & mesh.Sf())
-   //       + rhorAUf*fvc::ddtCorr(rho, U, phi)
+            fvc::flux(rhoeps*HbyA)
+    //      + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
-
-    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
     
-    // shamelessly stolen from interFoam
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
-    (
-        p.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-	  + phiUs.boundaryField()
-          - fvOptions.relative((mesh.Sf().boundaryField() & U.boundaryField()))
-	    * rhoeps.boundaryField()
-        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
-    );
+    // flux without pressure gradient contribution
+    phi = phiHbyA + phiUs;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rhoeps, U, phi, rhorAUf);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -86,48 +45,26 @@ else
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi*voidfraction, p)
-          + fvc::div(phiHbyA)
+          + fvc::div(phi)
           - fvm::laplacian(rhorAUf, p)
-	  + fvc::div(phiUs)
           ==
             fvOptions(psi, p, rho.name())
         );
 
-        fvOptions.constrain(pEqn);
-
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA + pEqn.flux()+phiUs;
-
-            // Explicitly relax pressure for momentum corrector
-            p.relax();
-
-            if (UfromPhi)
-	        U = fvc::reconstruct( phi/fvc::interpolate(rhoeps) );
-	    else
-                U = HbyA + rAU*(voidfraction*fvc::reconstruct((pEqn.flux())/rhorAUf) + Ksl*Us);
-	    
-            U.correctBoundaryConditions();
-            fvOptions.correct(U);
-	    
-divphi = fvc::div(phi);
-divphi2 = fvc::div(linearInterpolate(rho*U*voidfraction) & mesh.Sf());
-gP = fvc::grad(p);
-gP2 = fvc::reconstruct((pEqn.flux())/rhorAUf);
+            phi += pEqn.flux();
         }
     }
 }
 
-
-
 #include "rhoEqn.H"
-
 #include "compressibleContinuityErrsPU.H"
 
 // Explicitly relax pressure for momentum corrector
-//p.relax();
+p.relax();
 
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
@@ -137,10 +74,16 @@ rho.relax();
 Info<< "rho max/min : " << max(rho).value()
     << " " << min(rho).value() << endl;
 
-    // old version to get U from p; can have problems at boundaries where grad(p) != snGrad(p)
-//U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);
-//U.correctBoundaryConditions();
-//fvOptions.correct(U);
+if (modelType=="A")
+{
+    U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);
+}
+else
+{
+    U = HbyA - rAU*(fvc::grad(p)-Ksl*Us);
+}
+U.correctBoundaryConditions();
+fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
 if (thermo.dpdt())

src/finiteVolume/Make/files

@@ -3,7 +3,6 @@ derivedFvPatchFields = $(fvPatchFields)/derived
 
 general = cfdTools/general
 
-$(derivedFvPatchFields)/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchFields.C
 $(derivedFvPatchFields)/uniformFixedValueTube/uniformFixedValueTubeFvPatchFields.C
 
 LIB = $(CFDEM_LIB_DIR)/libfiniteVolumeCFDEM

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C

@@ -88,7 +88,7 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(DataEntry<Type>::New("uniformValue", dict)),
+    uniformValue_(Function1<Type>::New("uniformValue", dict)),
 	pName_("p"), //JOKER
 	phiName_("phi"), //JOKER
     velocityFieldName_("U"),

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H

@@ -40,7 +40,7 @@ SourceFiles
 
 #include "Random.H"
 #include "fixedValueFvPatchFields.H"
-#include "DataEntry.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ class uniformFixedValueTubeFvPatchField
 {
     // Private data
 
-        autoPtr<DataEntry<Type> > uniformValue_;
+        autoPtr<Function1<Type> > uniformValue_;
 
 		word pName_; //JOKER pressure
 		

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.C

@@ -73,7 +73,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(DataEntry<Type>::New("uniformValue", dict)),
+    uniformValue_(Function1<Type>::New("uniformValue", dict)),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.H

@@ -40,7 +40,7 @@ SourceFiles
 
 #include "Random.H"
 #include "fixedValueFvPatchFields.H"
-#include "DataEntry.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ class uniformFixedValueVoidfractionFvPatchField
 {
     // Private data
 
-        autoPtr<DataEntry<Type> > uniformValue_;
+        autoPtr<Function1<Type> > uniformValue_;
 
         word voidfractionFieldName_;
 

src/lagrangian/cfdemParticle/Make/files

@@ -118,7 +118,6 @@ $(locateModels)/turboEngineSearch/turboEngineSearch.C
 $(locateModels)/engineSearchMany2Many/engineSearchMany2Many.C
 $(locateModels)/engineSearchIB/engineSearchIB.C
 
-
 $(meshMotionModels)/meshMotionModel/meshMotionModel.C
 $(meshMotionModels)/meshMotionModel/newMeshMotionModel.C
 $(meshMotionModels)/noMeshMotion/noMeshMotion.C
@@ -139,7 +138,6 @@ $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
 $(dataExchangeModels)/twoWayFiles/twoWayFiles.C
 $(dataExchangeModels)/noDataExchange/noDataExchange.C
 $(dataExchangeModels)/twoWayMPI/twoWayMPI.C
-$(dataExchangeModels)/twoWayMany2Many/twoWayMany2Many.C
 
 $(averagingModels)/averagingModel/averagingModel.C
 $(averagingModels)/averagingModel/newAveragingModel.C

src/lagrangian/cfdemParticle/Make/options

@@ -12,7 +12,7 @@ EXE_INC = \
     -I ./cfdemParticle \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/OpenFOAM/containers/HashTables/labelHashSet \
@@ -25,8 +25,7 @@ LIB_LIBS = \
      $(PLIBS) \
     -L$(CFDEM_LIB_DIR) \
     -lfiniteVolume \
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lincompressibleTurbulenceModels \
     -lmeshTools \
     -llagrangian \
     -lmpi_cxx \

src/lagrangian/cfdemParticle/cfdTools/continuityErrorPhiPU.H

@@ -34,7 +34,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    volScalarField contErr( fvc::div(phiGes) + fvc::ddt(voidfraction) );
+    volScalarField contErr( fvc::div(phi) + fvc::ddt(voidfraction) );
 
     scalar sumLocalContErr = runTime.deltaTValue()*
         mag(contErr)().weightedAverage(mesh.V()).value();

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -127,15 +127,7 @@ cfdemCloud::cfdemCloud
     ),
     turbulence_
     (
-        #if defined(version21) || defined(version16ext)
-            #ifdef compre
-                mesh.lookupObject<compressible::turbulenceModel>
-            #else
-                mesh.lookupObject<incompressible::turbulenceModel>
-            #endif
-        #elif defined(version15)
-            mesh.lookupObject<incompressible::RASModel>
-        #endif
+        mesh.lookupObject<turbulenceModel>
         (
             turbulenceModelType_
         )
@@ -753,8 +745,8 @@ void cfdemCloud::resetArray(double**& array,int length,int width,double resetVal
 
 void cfdemCloud::otherForces(volVectorField& forcefield)
 {
-  forcefield.internalField() = vector::zero;
-  forcefield.boundaryField() = vector::zero;
+  forcefield.primitiveFieldRef() = vector::zero;
+  forcefield.boundaryFieldRef() = vector::zero;
   for (int i=0;i<otherForceModels_.size();i++)
       forcefield += otherForceModel_[i]().exportForceField();
 }

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -49,11 +49,7 @@ SourceFiles
 #include "fvCFD.H"
 #include "IFstream.H"
 
-#if defined(version21) || defined(version16ext)
-    #include <turbulenceModel.H>
-#elif defined(version15)
-    #include <RASModel.H>
-#endif
+#include <turbulenceModel.H>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -160,15 +156,7 @@ protected:
 
     mutable volScalarField ddtVoidfraction_;
 
-    #if defined(version21) || defined(version16ext)
-        #ifdef compre
-            const compressible::turbulenceModel& turbulence_;
-        #else
-            const incompressible::turbulenceModel& turbulence_;
-        #endif
-    #elif defined(version15)
-        const incompressible::RASModel& turbulence_;
-    #endif
+    const turbulenceModel& turbulence_;
 
     autoPtr<forceModel>* forceModel_;
 
@@ -371,16 +359,8 @@ public:
 
        inline autoPtr<liggghtsCommandModel>* liggghtsCommand() const;
 
-       #if defined(version21) || defined(version16ext)
-            #ifdef compre
-                inline const compressible::turbulenceModel& turbulence() const;
-            #else
-                inline const incompressible::turbulenceModel& turbulence() const;
-            #endif
-       #elif defined(version15)
-           inline const incompressible::RASModel& turbulence() const;
-       #endif
-
+       inline const turbulenceModel& turbulence() const;
+       
     // Write
 
       // write cfdemCloud internal data

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -302,15 +302,7 @@ inline autoPtr<liggghtsCommandModel>* cfdemCloud::liggghtsCommand() const
     return liggghtsCommand_;
 }
 
-#if defined(version21) || defined(version16ext)
-    #ifdef compre
-        inline const compressible::turbulenceModel& cfdemCloud::turbulence() const
-    #else
-        inline const incompressible::turbulenceModel& cfdemCloud::turbulence() const
-    #endif
-#elif defined(version15)
-    inline const incompressible::RASModel& cfdemCloud::turbulence() const
-#endif
+inline const turbulenceModel& cfdemCloud::turbulence() const
 {
     return turbulence_;
 }

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C

@@ -119,16 +119,16 @@ const thermCondModel& cfdemCloudEnergy::thermCondM()
 
 void cfdemCloudEnergy::energyContributions(volScalarField& Qsource)
 {
-    Qsource.internalField()=0.0;
-    Qsource.boundaryField()=0.0;
+    Qsource.primitiveFieldRef()=0.0;
+    Qsource.boundaryFieldRef()=0.0;
     for (int i=0;i<nrEnergyModels();i++)
         energyM(i).addEnergyContribution(Qsource);
 }
 
 void cfdemCloudEnergy::energyCoefficients(volScalarField& Qcoeff)
 {
-    Qcoeff.internalField()=0.0;
-    Qcoeff.boundaryField()=0.0;
+    Qcoeff.primitiveFieldRef()=0.0;
+    Qcoeff.boundaryFieldRef()=0.0;
     for (int i=0;i<nrEnergyModels();i++)
         energyM(i).addEnergyCoefficient(Qcoeff);
 }

src/lagrangian/cfdemParticle/etc/OFversion/OFversion.H

@@ -1,4 +1,5 @@
-#define version23  // currently being used
+#define version4x
+//#define version23  // currently being used
 //#define version22
 //#define version21
 //#define version16ext

src/lagrangian/cfdemParticle/etc/functions.sh

@@ -76,17 +76,17 @@ compileLib()
         i=$((${#str}-4))
         ending=${str:$i:4}
         if [[ $ending == "Comp" ]]; then
-                echo "Compiling a compressible library - so doing an rmdepall of incomp library first."
+                echo "Compiling a compressible library - so doing an wrmdep -a of incomp library first."
                 echo "Please make sure to have the compressible libraries first in the library-list.txt!"
                 cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle
                 echo "changing to $PWD"
-                rmdepall 2>&1 | tee -a $logpath/$logfileName
+                wrmdep -a 2>&1 | tee -a $logpath/$logfileName
                 cd $casePath
                 echo "changing to $PWD"
             else
                 echo "Compiling a incompressible library."
         fi
-        rmdepall 2>&1 | tee -a $logpath/$logfileName
+        wrmdep -a 2>&1 | tee -a $logpath/$logfileName
         wclean 2>&1 | tee -a $logpath/$logfileName
     #fi
     wmake libso 2>&1 | tee -a $logpath/$logfileName
@@ -128,7 +128,7 @@ compileSolver()
 
     #- wclean and wmake
     #if [ $doClean != "noClean" ]; then
-        rmdepall 2>&1 | tee -a $logpath/$logfileName
+        wrmdep -a 2>&1 | tee -a $logpath/$logfileName
         wclean 2>&1 | tee -a $logpath/$logfileName
     #fi
     
@@ -284,7 +284,7 @@ cleanCFDEM()
 
             cd  $path
             echo "cleaning library $PWD"
-            rmdepall
+            wrmdep -a
             wclean    
             rm -r ./Make/linux*
             rm -r ./lnInclude
@@ -338,7 +338,7 @@ cleanCFDEM()
 
             cd  $path            
             echo "cleaning solver $PWD"
-            rmdepall
+            wrmdep -a
             wclean    
     done
 }

src/lagrangian/cfdemParticle/etc/solver-list.txt

@@ -1,7 +1,5 @@
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
-cfdemSolverPisoTemp/dir
 cfdemSolverRhoPimple/dir
 cfdemSolverRhoPimpleChem/dir
-cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir

src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C

@@ -59,7 +59,7 @@ void averagingModel::undoVectorAverage
 {
 // WARNING - not sure if this is valid for dilute model!!!
 
-    if(!single) fieldPrev.internalField() = fieldNext.internalField();
+    if(!single) fieldPrev.ref() = fieldNext.ref();
 
     label cellI;
     vector valueVec;
@@ -303,13 +303,13 @@ void averagingModel::setDSauter
 
 void averagingModel::resetVectorAverage(volVectorField& prev,volVectorField& next,bool single) const
 {
-    if(!single) prev.internalField() = next.internalField();
-    next.internalField() = vector::zero;
+    if(!single) prev.ref() = next.ref();
+    next.primitiveFieldRef() = vector::zero;
 }
 
 void averagingModel::resetWeightFields() const
 {
-    UsWeightField_.internalField() = 0;
+    UsWeightField_.ref() = 0;
 }
 
 
@@ -352,12 +352,12 @@ tmp<volVectorField> Foam::averagingModel::UsInterp() const
 
     if (particleCloud_.dataExchangeM().couplingStep() > 1)
     {
-        tsource() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
+        tsource.ref() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
                     + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_;
     }
     else
     {
-        tsource() = UsNext_;
+        tsource.ref() = UsNext_;
     }
 
     return tsource;

src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C

@@ -246,13 +246,13 @@ void species::execute()
 
 
         // pull changeOfSpeciesMass_, transform onto fields changeOfSpeciesMassFields_, add them up on changeOfGasMassField_
-        changeOfGasMassField_.internalField() = 0.0;
-        changeOfGasMassField_.boundaryField() = 0.0;
+        changeOfGasMassField_.primitiveFieldRef() = 0.0;
+        changeOfGasMassField_.boundaryFieldRef() = 0.0;
         for (int i=0; i<speciesNames_.size();i++)
         {
             particleCloud_.dataExchangeM().getData(mod_spec_names_[i],"scalar-atom",changeOfSpeciesMass_[i]);
-            changeOfSpeciesMassFields_[i].internalField() = 0.0;
-            changeOfSpeciesMassFields_[i].boundaryField() = 0.0;
+            changeOfSpeciesMassFields_[i].primitiveFieldRef() = 0.0;
+            changeOfSpeciesMassFields_[i].boundaryFieldRef() = 0.0;
             particleCloud_.averagingM().setScalarSum
             (
                 changeOfSpeciesMassFields_[i],
@@ -262,7 +262,7 @@ void species::execute()
             );
 
             // take care for implementation in LIGGGHTS: species produced from particles defined positive
-            changeOfSpeciesMassFields_[i].internalField() /= changeOfSpeciesMassFields_[i].mesh().V();
+            changeOfSpeciesMassFields_[i].primitiveFieldRef() /= changeOfSpeciesMassFields_[i].mesh().V();
             changeOfSpeciesMassFields_[i].correctBoundaryConditions();
             changeOfGasMassField_ += changeOfSpeciesMassFields_[i];
             Info << "total conversion of species" << speciesNames_[i] << " = " << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V()) << endl;

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C

@@ -142,14 +142,14 @@ void heatTransferGunn::calcEnergyContribution()
     allocateMyArrays();
 
     // reset Scalar field
-    QPartFluid_.internalField() = 0.0;
+    QPartFluid_.primitiveFieldRef() = 0.0;
 
     // get DEM data
     particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_);
     
     if(calcPartTempField_)
     {       
-	partTempField_.internalField() = 0.0;
+	partTempField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
         particleCloud_.averagingM().setScalarAverage
         (
@@ -252,7 +252,7 @@ void heatTransferGunn::calcEnergyContribution()
         NULL
     );
 
-    QPartFluid_.internalField() /= -QPartFluid_.mesh().V();
+    QPartFluid_.primitiveFieldRef() /= -QPartFluid_.mesh().V();
 
     // limit source term
     forAll(QPartFluid_,cellI)

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnImplicit/heatTransferGunnImplicit.C

@@ -92,7 +92,7 @@ void heatTransferGunnImplicit::calcEnergyContribution()
         NULL
     );
 
-    QPartFluidCoeff_.internalField() /= -QPartFluidCoeff_.mesh().V();
+    QPartFluidCoeff_.primitiveFieldRef() /= -QPartFluidCoeff_.mesh().V();
     
     QPartFluidCoeff_.correctBoundaryConditions();
 

src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C

@@ -64,11 +64,7 @@ Archimedes::Archimedes
     propsDict_(dict.subDict(typeName + "Props")),
     twoDimensional_(false),
     gravityFieldName_(propsDict_.lookup("gravityFieldName")),
-    #if defined(version21) || defined(version16ext)
-        g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
-    #elif defined(version15)
-        g_(dimensionedVector(sm.mesh().lookupObject<IOdictionary>("environmentalProperties").lookup(gravityFieldName_)).value())
-    #endif
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
 {
 
     //Append the field names to be probed

src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.H

@@ -65,11 +65,7 @@ private:
 
     word gravityFieldName_;
 
-    #ifdef version21
-        const uniformDimensionedVectorField& g_;    // ref to gravity
-    #elif defined(version16ext) || defined(version15)
-        const dimensionedVector& g_;    // ref to gravity
-    #endif
+    const uniformDimensionedVectorField& g_;    // ref to gravity
 
 
 public:

src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C

@@ -66,11 +66,7 @@ ArchimedesIB::ArchimedesIB
     voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), //mod by alice
     voidfractions_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),//mod by alice
     gravityFieldName_(propsDict_.lookup("gravityFieldName")),
-    #if defined(version21) || defined(version16ext)
-        g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
-    #elif  defined(version15)
-        g_(dimensionedVector(sm.mesh().lookupObject<IOdictionary>("environmentalProperties").lookup(gravityFieldName_)).value())
-    #endif
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, "archimedesIBF.logDat");

src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.H

@@ -71,12 +71,8 @@ private:
 
     word gravityFieldName_;
 
-    #ifdef version21
-        const uniformDimensionedVectorField& g_;    // ref to gravity
-    #elif defined(version16ext) || defined(version15)
-        const dimensionedVector& g_;    // ref to gravity
-    #endif
-
+    const uniformDimensionedVectorField& g_;    // ref to gravity
+    
 
 public:
 

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C

@@ -304,7 +304,7 @@ void FinesFields::update()
 
 void FinesFields::calcSource() 
 {
-    Sds_.internalField()=0;
+    Sds_.primitiveFieldRef()=0;
     scalar f(0.0);
     scalar critpore(0.0);
     scalar deltaAlpha(0.0);
@@ -453,7 +453,7 @@ void FinesFields::updateDragCoeff()
 	        Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
 	    else
 	        Cd = 0.44;
-	    DragCoeff_.boundaryField()[patchI][faceI] = Cd * beta.boundaryField()[patchI][faceI];
+	    DragCoeff_.boundaryFieldRef()[patchI][faceI] = Cd * beta.boundaryFieldRef()[patchI][faceI];
         }
         
     DragCoeff_ = max( DragCoeff_, dimensionedScalar("SMALL", dimensionSet(1,-3,-1,0,0), SMALL) );

src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C

@@ -130,7 +130,7 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
     allocateMyArrays();
 
     // reset Scalar field
-    EuField.internalField() = 0.0;
+    EuField.primitiveFieldRef() = 0.0;
 
     // get DEM data
     particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_);
@@ -232,7 +232,7 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
     );
 
     // scale with -1/(Vcell*rho*Cp)
-    EuField.internalField() /= -rhoField.internalField()*Cp_*EuField.mesh().V();
+    EuField.primitiveFieldRef() /= -rhoField.internalField()*Cp_*EuField.mesh().V();
 
     // limit source term
     forAll(EuField,cellI)

src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C

@@ -138,8 +138,8 @@ void dSauter::setForce() const
 	}
     }
     
-    d2Field_.internalField() = 0.0;
-    d3Field_.internalField() = 0.0;
+    d2Field_.primitiveFieldRef() = 0.0;
+    d3Field_.primitiveFieldRef() = 0.0;
     
     particleCloud_.averagingM().setScalarSum
     (

src/lagrangian/cfdemParticle/subModels/forceModel/particleCellVolume/particleCellVolume.C

@@ -120,7 +120,7 @@ void particleCellVolume::setForce() const
     {
         if(verbose_) Info << "particleCellVolume.C - setForce()" << endl;
 
-        scalarField_.internalField()=0.;
+        scalarField_.ref()=0.;
 
         // get reference to actual field
         const volScalarField& field = mesh_.lookupObject<volScalarField>(scalarFieldName_);
@@ -143,8 +143,8 @@ void particleCellVolume::setForce() const
                 scalarField2_[cellI] = 0.;
             }
         }
-        scalarField_.internalField() = gSum(scalarField_);
-        scalarField2_.internalField() = gSum(scalarField2_);
+        scalarField_.ref() = gSum(scalarField_);
+        scalarField2_.ref() = gSum(scalarField2_);
 
         if(verbose_)
         {

src/lagrangian/cfdemParticle/subModels/forceModel/volWeightedAverage/volWeightedAverage.C

@@ -179,7 +179,7 @@ void volWeightedAverage::setForce() const
             MPI_Allreduce(&totVol, &totVol_all, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
             integralValue = gSum(scalarFields_[i]);
             volWeightedAverage = integralValue / (totVol_all+SMALL);
-            scalarFields_[i].internalField() = volWeightedAverage;
+            scalarFields_[i].ref() = volWeightedAverage;
 
             if(verbose_)
             {
@@ -224,7 +224,7 @@ void volWeightedAverage::setForce() const
 
             MPI_Allreduce(&totVol, &totVol_all, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
             volWeightedAverage = gSum(vectorFields_[i]) / (totVol_all+SMALL);
-            vectorFields_[i].internalField() = volWeightedAverage;
+            vectorFields_[i].ref() = volWeightedAverage;
 
             if(verbose_)
             {

src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.C

@@ -61,14 +61,8 @@ engineSearch::engineSearch
 :
     locateModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
-    //faceDecomp_(propsDict_.lookup("faceDecomp")),
     treeSearch_(propsDict_.lookup("treeSearch")),
-    #ifdef version16ext
-        searchEngine_(particleCloud_.mesh(),false) //(particleCloud_.mesh(),faceDecomp_)
-    #elif defined(version21)
-        searchEngine_(particleCloud_.mesh(),polyMesh::FACEPLANES) // FACEPLANES or FACECENTRETETS; FACEDIAGTETS not stable
-    #endif
-    //searchEngine_(particleCloud_.mesh(),faceDecomp_) // only 2.0.x
+    searchEngine_(particleCloud_.mesh(), polyMesh::FACE_PLANES)
 {}
 
 

src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.H

@@ -61,8 +61,6 @@ private:
 
     dictionary propsDict_;
 
-    //Switch faceDecomp_;
-
     Switch treeSearch_;
 
 protected:

src/lagrangian/cfdemParticle/subModels/locateModel/standardSearch/standardSearch.C

@@ -89,11 +89,7 @@ label standardSearch::findCell
             for(int i=0;i<3;i++) position[i] = positions[index][i];
 
             // find cell
-            #ifdef version16ext
-                cellIDs[index][0] = particleCloud_.mesh().findCell(position);
-            #elif defined(version21)
-                cellIDs[index][0] = particleCloud_.mesh().findCell(position, polyMesh::FACEPLANES);
-            #endif
+            cellIDs[index][0] = particleCloud_.mesh().findCell(position, polyMesh::FACE_PLANES);
         }
         else cellIDs[index][0] = -1;
     }
@@ -108,11 +104,7 @@ label standardSearch::findSingleCell
 ) const
 {
     // find cell
-    #ifdef version16ext
-        return particleCloud_.mesh().findCell(position);
-    #elif defined(version21)
-        return particleCloud_.mesh().findCell(position, polyMesh::FACEPLANES);
-    #endif
+    return particleCloud_.mesh().findCell(position, polyMesh::FACE_PLANES);
 }
 
 

src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearch/turboEngineSearch.C

@@ -64,11 +64,7 @@ turboEngineSearch::turboEngineSearch
     propsDict_(dict.subDict(typeName + "Props")),
     treeSearch_(propsDict_.lookup("treeSearch")),
     bb_(particleCloud_.mesh().points(),false),
-    #ifdef version16ext
-        searchEngine_(particleCloud_.mesh(),false) //(particleCloud_.mesh(),faceDecomp_)
-    #elif defined(version21)
-        searchEngine_(particleCloud_.mesh(),polyMesh::FACEPLANES) // FACEPLANES or FACECENTRETETS; FACEDIAGTETS not stable
-    #endif
+    searchEngine_(particleCloud_.mesh(), polyMesh::FACE_PLANES)
 {}
 
 

src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C

@@ -147,18 +147,18 @@ tmp<volVectorField> explicitCouple::expMomSource() const
                 if (magF > fLimit_[i]) fNext_[cellI][i] *= fLimit_[i]/magF;
             }
         }
-        tsource() = fPrev_;
+        tsource.ref() = fPrev_;
     }else
     {
-        tsource() = (1 - tsf) * fPrev_ + tsf * fNext_;
+        tsource.ref() = (1 - tsf) * fPrev_ + tsf * fNext_;
     }
     return tsource;
 }
 
 void Foam::explicitCouple::resetMomSourceField() const
 {
-    fPrev_.internalField() = fNext_.internalField();
-    fNext_.internalField() = vector::zero;
+    fPrev_.ref() = fNext_.ref();
+    fNext_.primitiveFieldRef() = vector::zero;
 }
 
 inline vector Foam::explicitCouple::arrayToField(label cellI) const

src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C

@@ -166,10 +166,10 @@ tmp<volScalarField> implicitCouple::impMomSource() const
             // limiter
             if (KslNext_[cellI] > KslLimit_) KslNext_[cellI] = KslLimit_;
         }
-        tsource() = KslPrev_;
+        tsource.ref() = KslPrev_;
     }else
     {
-        tsource() = (1 - tsf) * KslPrev_ + tsf * KslNext_;
+        tsource.ref() = (1 - tsf) * KslPrev_ + tsf * KslNext_;
     }
 
     return tsource;
@@ -177,8 +177,8 @@ tmp<volScalarField> implicitCouple::impMomSource() const
 
 void Foam::implicitCouple::resetMomSourceField() const
 {
-    KslPrev_.internalField() = KslNext_.internalField();
-    KslNext_.internalField() = 0;
+    KslPrev_.ref() = KslNext_.ref();
+    KslNext_.primitiveFieldRef() = 0;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/otherForceModel/expParticleForces/expParticleForces.C

@@ -74,7 +74,7 @@ tmp<volVectorField> expParticleForces::exportForceField()
         )
     );
     
-    volVectorField& source = tsource();
+    volVectorField& source = tsource.ref();
     
     // negative sign in sum because force on particles = - force on fluid
     for(int i=0; i<particleCloud_.nrMomCoupleModels(); i++)

src/lagrangian/cfdemParticle/subModels/otherForceModel/gravity/gravity.C

@@ -79,7 +79,7 @@ tmp<volVectorField> gravity::exportForceField()
         )
     );
     
-    volVectorField& source = tsource();
+    volVectorField& source = tsource.ref();
     
     source = rhoG_ * voidfraction_ * g_;
     

src/lagrangian/cfdemParticle/subModels/otherForceModel/weightSecondaryPhase/weightSecondaryPhase.C

@@ -83,7 +83,7 @@ tmp<volVectorField> weightSecondaryPhase::exportForceField()
         )
     );
     
-    volVectorField& source = tsource();
+    volVectorField& source = tsource.ref();
     
     source = rho_ * alpha_ * g_;
     

src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C

@@ -100,7 +100,7 @@ void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
     volScalarField sSmoothField = sSmoothField_;
 
     sSmoothField.dimensions().reset(fieldSrc.dimensions());
-    sSmoothField.internalField()=fieldSrc.internalField();
+    sSmoothField.ref()=fieldSrc.internalField();
     sSmoothField.correctBoundaryConditions();
     sSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     sSmoothField.oldTime()=fieldSrc;
@@ -141,7 +141,7 @@ void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
     volVectorField vSmoothField = vSmoothField_;
 
     vSmoothField.dimensions().reset(fieldSrc.dimensions());
-    vSmoothField.internalField()=fieldSrc.internalField();
+    vSmoothField.ref()=fieldSrc.internalField();
     vSmoothField.correctBoundaryConditions();
     vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     vSmoothField.oldTime()=fieldSrc;
@@ -176,7 +176,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
     volVectorField vSmoothField = vSmoothField_;
 
     vSmoothField.dimensions().reset(fieldSrc.dimensions());
-    vSmoothField.internalField()=fieldSrc.internalField();
+    vSmoothField.ref()=fieldSrc.internalField();
     vSmoothField.correctBoundaryConditions();
     vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     vSmoothField.oldTime()=fieldSrc;
@@ -210,7 +210,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
     forAll(vSmoothField,cellI)
     {
         if ( mag(vSmoothField.oldTime().internalField()[cellI]) > 0.0f)  // have a vector in the OLD vSmoothField, so keep it!
-            NLarge()[cellI] = sourceStrength;
+            NLarge.ref()[cellI] = sourceStrength;
     }
 
     // do the smoothing

src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.C

@@ -103,14 +103,14 @@ tmp<volScalarField> SyamlalThermCond::thermCond() const
 	    dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
         )
     );
-    volScalarField& svf = tvf();
+    volScalarField& svf = tvf.ref();
     
     svf = (1-sqrt(1-voidfraction_)) / (voidfraction_) * kf0_;
     
     // if a wallQFactor field is present, use it to scale heat transport through a patch
     if (hasWallQFactor_)
         forAll(wallQFactor_.boundaryField(), patchi)
-            svf.boundaryField()[patchi] *= wallQFactor_.boundaryField()[patchi];
+            svf.boundaryFieldRef()[patchi] *= wallQFactor_.boundaryField()[patchi];
   
     return tvf;
 }

src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C

@@ -91,7 +91,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
 {
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     scalar radius(-1);
     scalar volume(0);

src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C

@@ -96,7 +96,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {

src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C

@@ -90,7 +90,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
 {
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     scalar radius(-1);
     scalar volume(0);

src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C

@@ -127,19 +127,19 @@ tmp<volScalarField> Foam::voidFractionModel::voidFractionInterp() const
     scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
     if(1-tsf < 1e-4 && particleCloud_.dataExchangeM().couplingStep() > 1) //tsf==1
     {
-        tsource() = voidfractionPrev_;
+        tsource.ref() = voidfractionPrev_;
     }
     else
     {
-        tsource() = (1 - tsf) * voidfractionPrev_ + tsf * voidfractionNext_;
+        tsource.ref() = (1 - tsf) * voidfractionPrev_ + tsf * voidfractionNext_;
     }
     return tsource;
 }
 
 void Foam::voidFractionModel::resetVoidFractions() const
 {
-    voidfractionPrev_.internalField() = voidfractionNext_.internalField();
-    voidfractionNext_.internalField() = 1;
+    voidfractionPrev_.ref() = voidfractionNext_.ref();
+    voidfractionNext_.ref() = 1;
 }
 
 /*void Foam::voidFractionModel::undoVoidFractions(double**const& mask) const

src/lagrangian/cfdemParticleComp/Make/options

@@ -15,7 +15,7 @@ EXE_INC = \
     -I ../cfdemParticle/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
@@ -28,8 +28,7 @@ LIB_LIBS = \
      $(PLIBS) \
     -L$(CFDEM_LIB_DIR) \
     -lfiniteVolume \
-    -lcompressibleRASModels \
-    -lcompressibleLESModels \
+    -lcompressibleTurbulenceModels \
     -lfluidThermophysicalModels \
     -lmeshTools \
     -llagrangian \

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U

@@ -30,29 +30,7 @@ boundaryField
 
     inlet
     {
-        /*type            flowRateInletVelocity;
-        flowRate        0.001;
-        value           uniform (0 0 0);*/
-
-        /*type            fixedValue;
-        value           uniform (0 0 0.0001);*/
-
-        //type            zeroGradient;
-
-/*        type            groovyBC;
-        variables       "Uend=vector(0,0,0.02);tEnd=0.1;";
-        valueExpression "((time() < tEnd) ? Uend/tEnd*time():Uend)";
-        value           uniform (0 0 0);*/
-
-        /*// 2.0.x, 1.6,ext
-        type            timeVaryingUniformFixedValue;
-        fileName        "steps_0p1s";
-        outOfBounds     clamp;
-        value           uniform (0 0 0);*/
-
-        // 2.1.x
-		//type              uniformFixedValue;
-        type            uniformFixedValueVoidfraction;
+        type            uniformFixedValue;
         uniformValue     table
         (
             (0.000 (0 0 0.002))
@@ -62,9 +40,6 @@ boundaryField
 
     outlet
     {
-        /*type            fluxCorrectedVelocity; //inletOutlet;
-        value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);*/
         type            zeroGradient;
     }
 }

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;
-
-turbulence      off;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties

@@ -70,7 +70,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType "RASProperties";//"LESProperties";//
+turbulenceModelType "turbulenceProperties";
 
 //===========================================================================//
 // sub-model properties

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/steps_0p1s

@@ -1,22 +0,0 @@
-(
-    (0.000 (0 0 0.002))
-    (0.010 (0 0 0.002))
-    (0.011 (0 0 0.004))
-    (0.020 (0 0 0.004))
-    (0.021 (0 0 0.006))
-    (0.030 (0 0 0.006))
-    (0.031 (0 0 0.008))
-    (0.040 (0 0 0.008))
-    (0.041 (0 0 0.010))
-    (0.050 (0 0 0.010))
-    (0.051 (0 0 0.012))
-    (0.060 (0 0 0.012))
-    (0.061 (0 0 0.014))
-    (0.070 (0 0 0.014))
-    (0.071 (0 0 0.016))
-    (0.080 (0 0 0.016))
-    (0.081 (0 0 0.018))
-    (0.090 (0 0 0.018))
-    (0.091 (0 0 0.020))
-    (0.100 (0 0 0.020))
-)

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict

@@ -53,8 +53,7 @@ maxCo 0.1;
 libs ( "libfiniteVolumeCFDEM.so" );
 
 functions
-(
-
+{
     probes
     {
         type        probes;
@@ -94,76 +93,5 @@ functions
         outputControl   timeStep;//outputTime;
         outputInterval  1;
     }
-
-/*
-    probes_inletFlux
-    {
-        type            faceSource;
-        functionObjectLibs ("libfieldFunctionObjects.so");
-        enabled         true;
-        outputControl   timeStep;   //outputTime;
-        log             true;       // log to screen?
-        valueOutput     true;       // Write values at run-time output times?
-        source          patch;   // Type of face source:
-        sourceName      inlet;         // faceZone name, see below
-        operation       sum;
-        fields
-        (
-            phi
-        );
-    }
-
-    probes_outletFlux
-    {
-        type            faceSource;
-        functionObjectLibs ("libfieldFunctionObjects.so");
-        enabled         true;
-        outputControl   timeStep;   //outputTime;
-        log             true;       // log to screen?
-        valueOutput     true;       // Write values at run-time output times?
-        source          patch;   // Type of face source:
-        sourceName      outlet;         // faceZone name, see below
-        operation       sum;
-        fields
-        (
-            phi
-        );
-    }
-
-    probes_wallFlux
-    {
-        type            faceSource;
-        functionObjectLibs ("libfieldFunctionObjects.so");
-        enabled         true;
-        outputControl   timeStep;   //outputTime;
-        log             true;       // log to screen?
-        valueOutput     true;       // Write values at run-time output times?
-        source          patch;   // Type of face source:
-        sourceName      wall;         // faceZone name, see below
-        operation       sum;
-        fields
-        (
-            phi
-        );
-    }*/
-   /*pressureDrop
-   {
-     type patchAverage;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields
-     (
-         p
-     );
-     factor  1;
-   }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/U

@@ -30,40 +30,16 @@ boundaryField
 
     inlet
     {
-        /*type            flowRateInletVelocity;
-        flowRate        0.001;
-        value           uniform (0 0 0);*/
-
-        /*type            fixedValue;
-        value           uniform (0 0 0.0001);*/
-
-        //type            zeroGradient;
-
-/*        type            groovyBC;
-        variables       "Uend=vector(0,0,0.02);tEnd=0.1;";
-        valueExpression "((time() < tEnd) ? Uend/tEnd*time():Uend)";
-        value           uniform (0 0 0);*/
-
-        /*// 2.0.x, 1.6,ext
-        type            timeVaryingUniformFixedValue;
-        fileName        "steps_0p1s";
-        outOfBounds     clamp;
-        value           uniform (0 0 0);*/
-
-        // 2.1.x
         type            uniformFixedValue;
         uniformValue     table
         (
             (0.000 (0 0 0.002))
-            (0.100 (0 0 0.020))
+            (0.500 (0 0 0.020))
         );
     }
 
     outlet
     {
-        /*type            fluxCorrectedVelocity; //inletOutlet;
-        value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);*/
         type            zeroGradient;
     }
 }

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/p

@@ -16,7 +16,7 @@ FoamFile
 
 dimensions      [0 2 -2 0 0 0 0];
 
-internalField   uniform 1;
+internalField   uniform 0e0;
 
 boundaryField
 {
@@ -29,13 +29,9 @@ boundaryField
     inlet
     {
         type            zeroGradient;
-        //type            fixedValue;
-        //value           uniform 100000;
     }
     outlet
     {
-        //type            zeroGradient;
-
         type            fixedValue;
         value           $internalField;
     }

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;
-
-turbulence      off;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/couplingProperties

@@ -55,12 +55,12 @@ smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; //
 
 forceModels
 (
-    //GidaspowDrag
+    GidaspowDrag
     //BeetstraDrag
     //DiFeliceDrag
     gradPForce
     viscForce
-    KochHillDrag
+    //KochHillDrag
     //DEMbasedDrag
     //RongDrag
     //Archimedes
@@ -75,7 +75,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType "RASProperties";//"LESProperties";//
+turbulenceModelType "turbulenceProperties";//"LESProperties";//
 
 //===========================================================================//
 // sub-model properties
@@ -143,6 +143,7 @@ GidaspowDragProps
 {
     verbose true;
     velFieldName "U";
+    granVelFieldName "Us";
     voidfractionFieldName "voidfraction";
     interpolation true;
     phi 1;

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/octave/totalPressureDrop.m

@@ -9,7 +9,7 @@ rhoG = 10			% density in kg/m3
 %path = '../probes/0/p'; % 2.1.x
 path = '../postProcessing/probes/0/p'; % 2.2.x
 columns=22;
-headerlines=4;
+headerlines=23;%4;
 data = loaddata(path,columns,headerlines);
 data=transpose(data);
 [x,y]=size(data)

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/controlDict

@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.1;
+endTime         0.5;
 
 deltaT          0.0005;
 
@@ -49,10 +49,9 @@ adjustTimeStep  no;
 
 maxCo 0.1;
 
-//libs ( "libgroovyBC.so" "libfiniteVolumeCFDEM.so");
 
 functions
-(
+{
 
     probes
     {
@@ -94,24 +93,5 @@ functions
         outputInterval  1;
     }
 
-   /*pressureDrop
-   {
-     type patchAverage;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields
-     (
-         p
-     );
-     factor  1;
-   }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/0/U

@@ -23,34 +23,11 @@ boundaryField
 
     wall
     {
-        //type            fixedValue;
-        //value           uniform (0 0 0);
         type            slip;
     }
 
     inlet
     {
-        /*type            flowRateInletVelocity;
-        flowRate        0.001;
-        value           uniform (0 0 0);*/
-
-        /*type            fixedValue;
-        value           uniform (0 0 0.0001);*/
-
-        //type            zeroGradient;
-
-        /*type            groovyBC;
-        variables       "Uend=vector(0,0,2);tEnd=0.1;";
-        valueExpression "((time() < tEnd) ? Uend/tEnd*time():Uend)";
-        value           uniform (0 0 0);*/
-
-        /*// 2.0.x, 1.6,ext
-        type            timeVaryingUniformFixedValue;
-        fileName        "steps_0p1s";
-        outOfBounds     clamp;
-        value           uniform (0 0 0);*/
-
-        // 2.1.x
         type            uniformFixedValue;
         uniformValue     table
         (
@@ -79,9 +56,6 @@ boundaryField
 
     outlet
     {
-        /*type            fluxCorrectedVelocity; //inletOutlet;
-        value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);*/
         type            zeroGradient;
     }
 }

tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;
-
-turbulence      off;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/couplingProperties

@@ -67,7 +67,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType "RASProperties";//"LESProperties";//
+turbulenceModelType "turbulenceProperties";
 
 //===========================================================================//
 // sub-model properties

tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/system/controlDict

@@ -149,8 +149,7 @@ DimensionedConstants
 }
 
 functions
-(
-
+{
     probes
     {
         type        probes;
@@ -190,25 +189,5 @@ functions
         outputControl   timeStep;//outputTime;
         outputInterval  1;
     }
-
-   /*pressureDrop
-   {
-     type patchAverage;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields
-     (
-         p
-     );
-     factor  1;
-   }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/CFD/0/U

@@ -23,65 +23,21 @@ boundaryField
 
     wall
     {
-        //type            fixedValue;
-        //value           uniform (0 0 0);
         type            slip;
     }
 
     inlet
     {
-        /*type            flowRateInletVelocity;
-        flowRate        0.001;
-        value           uniform (0 0 0);*/
-
-        /*type            fixedValue;
-        value           uniform (0 0 0.0001);*/
-
-        //type            zeroGradient;
-
-/*        type            groovyBC;
-        variables       "Uend=vector(0,0,0.02);tEnd=0.1;";
-        valueExpression "((time() < tEnd) ? Uend/tEnd*time():Uend)";
-        value           uniform (0 0 0);*/
-
-        /*// 2.0.x, 1.6,ext
-        type            timeVaryingUniformFixedValue;
-        fileName        "steps_0p1s";
-        outOfBounds     clamp;
-        value           uniform (0 0 0);*/
-
-        // 2.1.x
         type            uniformFixedValue;
         uniformValue     table
         (
             (0.000 (0 0 0.002))
-            /*(0.010 (0 0 0.002))
-            (0.011 (0 0 0.004))
-            (0.020 (0 0 0.004))
-            (0.021 (0 0 0.006))
-            (0.030 (0 0 0.006))
-            (0.031 (0 0 0.008))
-            (0.040 (0 0 0.008))
-            (0.041 (0 0 0.010))
-            (0.050 (0 0 0.010))
-            (0.051 (0 0 0.012))
-            (0.060 (0 0 0.012))
-            (0.061 (0 0 0.014))
-            (0.070 (0 0 0.014))
-            (0.071 (0 0 0.016))
-            (0.080 (0 0 0.016))
-            (0.081 (0 0 0.018))
-            (0.090 (0 0 0.018))
-            (0.091 (0 0 0.020))*/
             (0.100 (0 0 0.020))
         );
     }
 
     outlet
     {
-        /*type            fluxCorrectedVelocity; //inletOutlet;
-        value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);*/
         type            zeroGradient;
     }
 }

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;
-
-turbulence      off;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/CFD/constant/couplingProperties_restart

@@ -69,7 +69,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType "RASProperties";//"LESProperties";//
+turbulenceModelType "turbulenceProperties";
 
 //===========================================================================//
 // sub-model properties

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/CFD/constant/couplingProperties_run

@@ -69,7 +69,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType "RASProperties";//"LESProperties";//
+turbulenceModelType "turbulenceProperties";
 
 //===========================================================================//
 // sub-model properties

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/CFD/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/CFD/system/controlDict

@@ -52,8 +52,7 @@ maxCo 0.1;
 //libs ( "libgroovyBC.so" "libfiniteVolumeCFDEM.so");
 
 functions
-(
-
+{
     probes
     {
         type        probes;
@@ -113,5 +112,5 @@ functions
      );
      factor  1;
    }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/CFD/system/controlDict_restart

@@ -52,8 +52,7 @@ maxCo 0.1;
 //libs ( "libgroovyBC.so" "libfiniteVolumeCFDEM.so");
 
 functions
-(
-
+{
     probes
     {
         type        probes;
@@ -113,5 +112,5 @@ functions
      );
      factor  1;
    }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/CFD/system/controlDict_run

@@ -52,8 +52,7 @@ maxCo 0.1;
 //libs ( "libgroovyBC.so" "libfiniteVolumeCFDEM.so");
 
 functions
-(
-
+{
     probes
     {
         type        probes;
@@ -113,5 +112,5 @@ functions
      );
      factor  1;
    }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPiso/settlingTestMPI/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;//kEpsilon;
-
-turbulence      on;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //

tutorials/cfdemSolverPiso/settlingTestMPI/CFD/constant/couplingProperties

@@ -69,7 +69,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType RASProperties;//LESProperties;//
+turbulenceModelType "turbulenceProperties";
 
 //===========================================================================//
 // sub-model properties

tutorials/cfdemSolverPiso/settlingTestMPI/CFD/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;//LESModel; //
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/cfdemSolverPisoMS/ErgunTestMPI/CFD/0/U

@@ -23,35 +23,11 @@ boundaryField
 
     wall
     {
-        //type            fixedValue;
-        //value           uniform (0 0 0);
         type            slip;
     }
 
     inlet
     {
-        /*type            flowRateInletVelocity;
-        flowRate        0.001;
-        value           uniform (0 0 0);*/
-
-/*        type            fixedValue;
-        value           uniform (0 0 0.1);*/
-
-        //type            zeroGradient;
-
-        /*// superficial velocity BC
-        type            groovyBC;
-        variables       "Usup=vector(0,0,2);alpha=voidfraction;tEnd=0.1;"; // should be used with zeroGradient voidfraction
-        valueExpression "((time() < tEnd) ? Usup/alpha*(time()/tEnd) :  Usup/alpha)";
-        value           uniform (0 0 0);*/
-
-        /*// 2.0.x, ext
-        type            timeVaryingUniformFixedValue;
-        fileName        "steps_0p1s";
-        outOfBounds     clamp;
-        value           uniform (0 0 0);*/
-
-        // 2.1.x
         type            uniformFixedValue;
         uniformValue     table
         (
@@ -62,9 +38,6 @@ boundaryField
     }
     outlet
     {
-        /*type            fluxCorrectedVelocity; //inletOutlet;
-        value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);*/
         type            zeroGradient;
     }
 }

tutorials/cfdemSolverPisoMS/ErgunTestMPI/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;
-
-turbulence      off;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //

tutorials/cfdemSolverPisoMS/ErgunTestMPI/CFD/constant/couplingProperties

@@ -70,7 +70,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType RASProperties;//LESProperties;//
+turbulenceModelType "turbulenceProperties";
 
 //===========================================================================//
 // sub-model properties

tutorials/cfdemSolverPisoMS/ErgunTestMPI/CFD/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/cfdemSolverPisoMS/ErgunTestMPI/CFD/system/controlDict

@@ -52,8 +52,7 @@ maxCo 0.1;
 //libs ("libOpenFOAM.so" "libgroovyBC.so");
 
 functions
-(
-
+{
     probes
     {
         type        probes;
@@ -93,25 +92,5 @@ functions
         outputControl   timeStep;//outputTime;
         outputInterval  1;
     }
-
-   /*pressureDrop
-   {
-     type patchAverage;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields
-     (
-         p
-     );
-     factor  1;
-   }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPisoScalar/packedBedTemp/CFD/constant/couplingProperties

@@ -62,7 +62,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType RASProperties;//LESProperties;//
+turbulenceModelType "turbulenceProperties";
 
 //===========================================================================//
 // sub-model properties

tutorials/cfdemSolverPisoScalar/packedBedTemp/CFD/system/controlDict

@@ -54,8 +54,7 @@ maxDeltaT       1;
 //libs ( "libgroovyBC.so" );
 
 functions
-(
-
+{
     probes
     {
         type        probes;
@@ -76,84 +75,5 @@ functions
         outputControl   timeStep;//outputTime;
         outputInterval  1;
     }
-
-/*
-// simpleFunctionObjects
-   heatFlux
-   {
-     type patchHeatFlux;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields    // name of temp field
-     (
-         T
-     );
-     cp 1007;       // cp in [J/(kg*K)]
-     factor  1.188; // density for incomp solvers!
-   }
-
-   massFlux
-   {
-     type patchMassFlow;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     factor  1.188; // density for incomp solvers!
-   }
-
-   pressureDrop
-   {
-     type patchAverage;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields
-     (
-         p
-     );
-     factor  1;
-   }
-
-   T
-   {
-     type patchAverage;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields    // name of temp field
-     (
-         T
-     );
-     factor  1;
-   }*/
-);
+};
 // ************************************************************************* //