Merge pull request #30 from ParticulateFlow/release

Release 17.08

.gitignore

@@ -5,4 +5,4 @@
 log_*
 log.*
 *~
-
+**/linux64GccDPInt32Opt

README

@@ -5,7 +5,8 @@
     www.cfdem.com
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
-                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright 2012-2015 DCS Computing GmbH, Linz
+                                Copyright 2015-     JKU Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -32,7 +33,7 @@ Description
 
 
 CFDEM(R) coupling provides an open source parallel coupled CFD-DEM framework 
-combining the strengths of LIGGGHTS(R) DEM code and the Open Source 
+combining the strengths of the LIGGGHTS(R) DEM code and the Open Source 
 CFD package OpenFOAM(R)(*). The CFDEM(R)coupling toolbox allows to expand 
 standard CFD solvers of OpenFOAM(R)(*) to include a coupling to the DEM 
 code LIGGGHTS(R). In this toolbox the particle representation within the 
@@ -46,8 +47,6 @@ Features are:
 
 - its modular approach allows users to easily implement new models
 - its MPI parallelization enables to use it for large scale problems
-- the "forum"_lws on CFD-DEM gives the possibility to exchange with other 
-  users / developers
 - the use of GIT allows to easily update to the latest version
 - basic documentation is provided
 
@@ -75,7 +74,7 @@ Core development of CFDEM(R)coupling is done by
 Christoph Goniva and Christoph Kloss, both at DCS Computing GmbH, 2012
 
 
-\*---------------------------------------------------------------------------*/
-(*) "OpenFOAM(R)"_of is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.
+/*---------------------------------------------------------------------------*\
+(*) "OpenFOAM(R)" is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.
 This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM® and OpenCFD® trade marks.
 \*---------------------------------------------------------------------------*/

applications/solvers/cfdemSolverIB/Make/options

@@ -1,24 +1,32 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -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 \
-    -l$(CFDEM_LIB_NAME)
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
 

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,11 +53,6 @@ Contributions
 
 #include "cellSet.H"
 
-#if defined(version22)
-    #include "meshToMeshNew.H"
-    #include "fvIOoptionList.H"
-#endif
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -67,14 +63,14 @@ int main(int argc, char *argv[])
 
     #include "createDynamicFvMesh.H"
 
+    #include "createControl.H"
+
+    #include "createTimeControls.H"
+
     #include "createFields.H"
 
     #include "initContinuityErrs.H"
 
-    #if defined(version22)
-    #include "createFvOptions.H"
-    #endif
-
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
     cfdemCloudIB particleCloud(mesh);
@@ -91,8 +87,9 @@ int main(int argc, char *argv[])
         interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
         mesh.update(); //dyM
 
-        #include "readPISOControls.H"
+        #include "readTimeControls.H"
         #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         // do particle stuff
         Info << "- evolve()" << endl;
@@ -107,45 +104,29 @@ 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.correctNonOrthogonal())
                 {
                     // Pressure corrector
 
@@ -156,20 +137,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();
                     }
@@ -188,10 +158,6 @@ int main(int argc, char *argv[])
         volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
         particleCloud.calcVelocityCorrection(p,U,phiIB,voidfractionNext);
 
-        #if defined(version22)
-        fvOptions.correct(U);
-        #endif
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/cfdemSolverPiso/Make/options

@@ -1,15 +1,23 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -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 \
-    -l$(CFDEM_LIB_NAME)
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPiso/UEqn.H

@@ -0,0 +1,23 @@
+fvVectorMatrix UEqn
+(
+    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
+  + particleCloud.divVoidfractionTau(U, voidfraction)
+  ==
+  - fvm::Sp(Ksl/rho,U)
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
+    fvOptions.correct(U);		    
+}
+else if (piso.momentumPredictor())
+{
+    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+    fvOptions.correct(U);
+}

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 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
-
-                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/cfdemSolverPiso/pEqn.H

@@ -0,0 +1,58 @@
+volScalarField rAU = 1.0/UEqn.A();
+
+surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
+
+volScalarField rAUvoidfraction("(voidfraction2|A(U))",rAU*voidfraction);
+
+if (modelType=="A")
+    rAUvoidfraction *= voidfraction;
+
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(voidfraction)*fvc::flux(HbyA) );
+
+volVectorField Uvoidfraction("(Uvoidfraction)", U * voidfraction);
+
+surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
+
+phi = phiHbyA + rAUf*(fvc::interpolate(Ksl/rho) * phiS);
+
+// rotating frames of references not tested yet
+MRF.makeRelative(phi);
+
+// adjustment of phi (only in cases w.o. p boundary conditions) not tested yet
+adjustPhi(phi, U, p);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, Uvoidfraction, phiHbyA, rAUvoidfraction, MRF);
+
+
+// Non-orthogonal pressure corrector loop
+while (piso.correctNonOrthogonal())
+{
+    // Pressure corrector
+  
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rAUvoidfraction, p) == fvc::div(phi) + particleCloud.ddtVoidfraction()
+    );
+    
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+    if (piso.finalNonOrthogonalIter())
+    {
+        phi -= pEqn.flux();
+    }
+}
+
+#include "continuityErrorPhiPU.H"
+
+if (modelType=="B" || modelType=="Bfull")
+    U = HbyA - rAU*fvc::grad(p) + Ksl/rho*Us*rAU;
+else
+    U = HbyA - voidfraction*rAU*fvc::grad(p) + Ksl/rho*Us*rAU;
+
+U.correctBoundaryConditions();
+fvOptions.correct(U);

applications/solvers/cfdemSolverPisoMS/Make/options

@@ -1,16 +1,24 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -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../cfdemSolverPiso \
+    -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 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
-
-            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,15 +1,24 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -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 \
-    -l$(CFDEM_LIB_NAME)
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoScalar/TEqn.H

@@ -0,0 +1,15 @@
+ // 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();

applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C

@@ -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 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
-
-                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

@@ -0,0 +1,59 @@
+// contributions to internal energy equation can be found in
+// Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998
+{
+    // dim he = J / kg
+    volScalarField& he = thermo.he();
+    particleCloud.energyContributions(Qsource);
+    particleCloud.energyCoefficients(QCoeff);
+
+    //thDiff=particleCloud.thermCondM().thermDiff();   
+    thCond=particleCloud.thermCondM().thermCond(); 
+
+    addSource = fvc::ddt(rhoeps, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        );
+
+    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
+
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rhoeps, he) + fvm::div(phi, he)
+      + addSource
+     // net heat transfer from particles to fluid
+      - Qsource
+      - fvm::Sp(QCoeff/Cpv, he)
+     // thermal conduction of the fluid with effective conductivity
+     // - fvm::laplacian(rhoeps*thDiff,he)
+      - fvm::laplacian(voidfraction*thCond/Cpv,he)
+      // + particle-fluid energy transfer due to work
+      // + fluid energy dissipation due to shearing
+     ==
+        fvOptions(rho, he)
+    );
+    
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+
+    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+    particleCloud.clockM().start(31,"postFlow");
+    particleCloud.postFlow();
+    particleCloud.clockM().stop("postFlow");
+}

applications/solvers/cfdemSolverRhoPimple/Make/files

@@ -0,0 +1,3 @@
+cfdemSolverRhoPimple.C
+
+EXE=$(CFDEM_APP_DIR)/cfdemSolverRhoPimple

applications/solvers/cfdemSolverRhoPimple/Make/options

@@ -0,0 +1,32 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+PFLAGS+= -Dcompre
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -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$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -l$(CFDEM_LIB_COMP_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverRhoPimple/UEqn.H

@@ -0,0 +1,33 @@
+// Solve the Momentum equation
+particleCloud.otherForces(fOther);
+
+tmp<fvVectorMatrix> tUEqn
+(
+    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();
+
+fvOptions.constrain(UEqn);
+
+if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    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);
+
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}

applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+Application
+    cfdemSolverRhoPimple
+
+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) 4.x,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "bound.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+
+#include "cfdemCloudEnergy.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+#include "thermCondModel.H"
+#include "energyModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudEnergy particleCloud(mesh);
+    #include "checkModelType.H"
+
+    turbulence->validate();
+  //        #include "compressibleCourantNo.H"
+  //  #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
+        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();
+        Ksl.correctBoundaryConditions();
+
+        //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");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        if (pimple.nCorrPIMPLE() <= 1)
+        {
+            #include "rhoEqn.H"
+        }
+
+        volScalarField rhoeps("rhoeps",rho*voidfraction);
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                // besides this pEqn, OF offers a "pimple consistent"-option
+                #include "pEqn.H"
+                rhoeps=rho*voidfraction;
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverRhoPimple/createFieldRefs.H

@@ -0,0 +1,2 @@
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();

applications/solvers/cfdemSolverRhoPimple/createFields.H

@@ -0,0 +1,253 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<psiThermo> pThermo
+    (
+        psiThermo::New(mesh)
+    );
+    psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
+    volScalarField& p = thermo.p();
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+           "addSource",
+           runTime.timeName(),
+           mesh,
+           IOobject::MUST_READ,
+           IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating fluid-particle heat flux field\n" << endl;
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl;
+    volScalarField QCoeff
+    (
+        IOobject
+        (
+            "QCoeff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    );
+
+  /*  Info<< "\nCreating thermal diffusivity field\n" << endl;
+    volScalarField thDiff
+    (
+        IOobject
+        (
+            "thDiff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
+    );
+  */
+    Info<< "\nCreating thermal conductivity field\n" << endl;
+    volScalarField thCond
+    (
+        IOobject
+        (
+            "thCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating heat capacity field\n" << endl;
+    volScalarField Cpv
+    (
+        IOobject
+        (
+            "Cpv",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
+    Info<< "Reading/calculating face flux field phi\n" << endl;
+    surfaceScalarField phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+         linearInterpolate(rho*U*voidfraction) & mesh.Sf()
+    );
+
+    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
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    mesh.setFluxRequired(p.name());
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "\nReading momentum exchange field Ksl\n" << endl;
+    volScalarField Ksl
+    (
+        IOobject
+        (
+            "Ksl",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+        //dimensionedScalar("0", dimensionSet(1, -3, -1, 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
+    );
+
+//===============================

applications/solvers/cfdemSolverRhoPimple/pEqn.H

@@ -0,0 +1,92 @@
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+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)
+{
+    tUEqn.clear();
+}
+
+if (pimple.transonic())
+{
+//  transonic version not implemented yet
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            fvc::flux(rhoeps*HbyA)
+    //      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+        )
+    );
+    
+    // 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())
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi*voidfraction, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrsPU.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
+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())
+{
+    dpdt = fvc::ddt(voidfraction,p);
+}

applications/solvers/cfdemSolverRhoPimple/rhoEqn.H

@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(voidfraction,rho)
+      + fvc::div(phi)
+      ==
+        fvOptions(rho)
+    );
+
+    fvOptions.constrain(rhoEqn);
+
+    rhoEqn.solve();
+
+    fvOptions.correct(rho);
+}
+
+// ************************************************************************* //

applications/utilities/cfdemPostproc/Make/options

@@ -1,5 +1,9 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -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 \
@@ -9,9 +13,11 @@ EXE_INC = \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
     -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
 

applications/utilities/cfdemPostproc/cfdemPostproc.C

@@ -33,7 +33,7 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
 
 #include "cfdemCloud.H"
 #include "dataExchangeModel.H"
@@ -72,8 +72,8 @@ int main(int argc, char *argv[])
     double **particleWeights_;
     double **particleVolumes_;
     double **particleV_;
-    double **cellIDs_;
-    
+    int **cellIDs_;
+
     particleCloud.dataExchangeM().allocateArray(positions_,0.,3);
     particleCloud.dataExchangeM().allocateArray(velocities_,0.,3);
     particleCloud.get_radii(radii_);  // get ref to radii
@@ -98,7 +98,7 @@ int main(int argc, char *argv[])
         particleCloud.averagingM().resetWeightFields();
         particleCloud.momCoupleM(0).resetMomSourceField();
 
-        particleCloud.dataExchangeM().couple();
+        particleCloud.dataExchangeM().couple(0);
 
         particleCloud.dataExchangeM().getData("x","vector-atom",positions_,count);
         particleCloud.dataExchangeM().getData("v","vector-atom",velocities_,count);
@@ -109,7 +109,7 @@ int main(int argc, char *argv[])
 
         particleCloud.voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_,particleV_);
 
-        voidfraction.internalField() = particleCloud.voidFractionM().voidFractionInterp();
+        voidfraction.ref() = particleCloud.voidFractionM().voidFractionInterp();
         voidfraction.correctBoundaryConditions();
 
         particleCloud.averagingM().setVectorAverage

applications/utilities/cfdemPostproc/createFields.H

@@ -10,7 +10,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("0", dimensionSet(0, 2, -2, 0, 0), 1.0)
+        dimensionedScalar("1", dimensionSet(0, 2, -2, 0, 0), 1.0)
     );
 
     Info<< "Reading physical velocity field U" << endl;
@@ -26,15 +26,29 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedVector("0", dimensionSet(0, 1, -1, 0, 0), vector::zero)
+        dimensionedVector("0", dimensionSet(0, 1, -1, 0, 0), Foam::vector::zero)
     );
 
 //========================
 // drag law modelling
 //========================
 
+    Info<< "Creating dummy density field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("1", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
 
-    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    Info<< "Reading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (
         IOobject
@@ -46,7 +60,7 @@
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("0", dimensionSet(0, 0, 0, 0, 0), 1.)
+        dimensionedScalar("1", dimensionSet(0, 0, 0, 0, 0), 1.0)
     );
 
 
@@ -62,7 +76,7 @@
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedVector("0", dimensionSet(0, 1, -1, 0, 0), vector::zero)
+        dimensionedVector("0", dimensionSet(0, 1, -1, 0, 0), Foam::vector::zero)
     );
 
 //========================

etc/OFversion/OFversion.H

@@ -0,0 +1,2 @@
+#define version4x
+

etc/additionalLibs

@@ -0,0 +1,8 @@
+# paths for additional libraries
+CFDEM_ADD_LIB_PATHS = \
+
+# additional libraries to be linked to solvers
+CFDEM_ADD_LIBS = \
+
+# additional static libraries to be linked to lagrangian library
+CFDEM_ADD_STATICLIBS = \

etc/bashrc

@@ -0,0 +1,259 @@
+#----------------------------------*-sh-*--------------------------------------
+# CFDEMcoupling
+#
+# Copyright 2009-2012 JKU Linz
+# Copyright 2012-2015 DCS Computing GmbH, Linz
+# Copyright 2015-     JKU Linz
+#------------------------------------------------------------------------------
+#
+# File
+#     etc/bashrc
+#
+# Description
+#     Startup file for CFDEMcoupling
+#     Sourced from ~/.profile or ~/.bashrc
+#
+# you can test the correctness using cfdemSystemTest.sh
+#------------------------------------------------------------------------------
+
+export CFDEM_PROJECT=CFDEM
+export CFDEM_VERSION=17.08
+
+################################################################################
+# USER EDITABLE PART: Changes made here may be lost with the next upgrade
+#
+# Please set to the appropriate path if the default is not correct.
+#
+# activate compatible OpenFOAM version
+. $HOME/OpenFOAM/OpenFOAM-4.x/etc/bashrc
+#
+# CFDEMcoupling
+export CFDEM_INST_DIR=$HOME/$CFDEM_PROJECT
+export CFDEM_PROJECT_DIR=$CFDEM_INST_DIR/CFDEMcoupling
+export CFDEM_PROJECT_USER_DIR=$CFDEM_INST_DIR/$USER-$WM_PROJECT_VERSION
+#
+# LIGGGHTS
+export CFDEM_LIGGGHTS_INST_DIR=$CFDEM_INST_DIR/LIGGGHTS
+#
+# LPP installation path
+export CFDEM_LPP_INST_DIR=$CFDEM_INST_DIR/LPP
+#
+# Path to additional libraries
+export CFDEM_ADD_LIBS_DIR=
+#
+# END OF (NORMAL) USER EDITABLE PART
+################################################################################
+
+#- The old dirs to be cleaned from the environment variables
+cfdemOldDirs="$CFDEM_LIGGGHTS_BIN_DIR $CFDEM_APP_DIR $CFDEM_LIB_DIR"
+
+export CFDEM_LIGGGHTS_SRC_DIR=$CFDEM_LIGGGHTS_INST_DIR/src
+export CFDEM_LIGGGHTS_BIN_DIR=$CFDEM_LIGGGHTS_INST_DIR/src-build
+
+#------------------------------------------------------------------------------
+# Source initialization functions
+#------------------------------------------------------------------------------
+. $WM_PROJECT_DIR/etc/config.sh/functions
+
+#------------------------------------------------------------------------------
+# Evaluate command-line parameters
+# these can be used to set/unset values
+#------------------------------------------------------------------------------
+_foamEval $@
+
+#------------------------------------------------------------------------------
+
+export CFDEM_SRC_DIR=$CFDEM_PROJECT_DIR/src
+export CFDEM_SOLVER_DIR=$CFDEM_PROJECT_DIR/applications/solvers
+export CFDEM_UT_DIR=$CFDEM_PROJECT_DIR/applications/utilities
+export CFDEM_DOC_DIR=$CFDEM_PROJECT_DIR/doc
+export CFDEM_TUT_DIR=$CFDEM_PROJECT_DIR/tutorials
+export CFDEM_LPP_DIR=$CFDEM_LPP_INST_DIR/src
+
+#------------------------------------------------------------------------------
+
+#- CFDEM lib name
+export CFDEM_LIB_NAME=lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#- CFDEM compressible lib name
+export CFDEM_LIB_COMP_NAME=lagrangianCFDEMcomp-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#- Check if additional libraries should be compiled together with solvers
+if [[ $CFDEM_ADD_LIBS_DIR == "" ]]; then
+    export CFDEM_ADD_LIBS_DIR=$CFDEM_PROJECT_DIR/etc
+else
+    echo "using CFDEM_ADD_LIBS_DIR=$CFDEM_ADD_LIBS_DIR defined by user."
+fi
+
+#------------------------------------------------------------------------------
+
+#- LMP Many2Many lib path and makefile
+export CFDEM_Many2ManyLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
+export CFDEM_Many2ManyLIB_MAKEFILENAME=fedora_fpic
+
+#------------------------------------------------------------------------------
+
+#- path to test harness
+export CFDEM_TEST_HARNESS_PATH=$CFDEM_PROJECT_USER_DIR/log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#- path to libraries
+export CFDEM_LIB_DIR=$CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/lib
+
+#- path to apps
+export CFDEM_APP_DIR=$CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/bin
+
+#- create directories
+mkdir -p $CFDEM_LIB_DIR
+mkdir -p $CFDEM_APP_DIR
+
+#- path to OF version flag file
+export CFDEM_OFVERSION_DIR=$CFDEM_PROJECT_DIR/etc/OFversion
+
+#------------------------------------------------------------------------------
+# Clean standard environment variables (PATH, LD_LIBRARY_PATH)
+#------------------------------------------------------------------------------
+foamClean=$WM_PROJECT_DIR/bin/foamCleanPath
+
+#- Clean PATH
+cleaned=`$foamClean "$PATH" "$cfdemOldDirs"` && PATH="$cleaned"
+
+#- Clean LD_LIBRARY_PATH
+cleaned=`$foamClean "$LD_LIBRARY_PATH" "$cfdemOldDirs"` \
+    && LD_LIBRARY_PATH="$cleaned"
+
+export PATH LD_LIBRARY_PATH
+
+#- add binary directories to $PATH
+_foamAddPath $CFDEM_APP_DIR:$CFDEM_LIGGGHTS_BIN_DIR
+_foamAddLib $CFDEM_LIB_DIR
+
+#------------------------------------------------------------------------------
+# Cleanup environment:
+#------------------------------------------------------------------------------
+unset cleaned foamClean cfdemOldDirs
+
+#------------------------------------------------------------------------------
+# Unload initialization functions:
+#------------------------------------------------------------------------------
+. $WM_PROJECT_DIR/etc/config.sh/functions
+
+#------------------------------------------------------------------------------
+# Settings for lpp postproc tool
+#------------------------------------------------------------------------------
+
+#- nr of procs for lpp tool
+export CFDEM_LPP_NPROCS=4
+
+#- nr of procs for lpp tool
+export CFDEM_LPP_CHUNKSIZE=1
+
+#- shortcut to run lpp
+alias lpp='python -i $CFDEM_LPP_DIR/lpp.py --cpunum $CFDEM_LPP_NPROCS --chunksize $CFDEM_LPP_CHUNKSIZE'
+
+#------------------------------------------------------------------------------
+# aliases for easy navigation (no changes necessary)
+#------------------------------------------------------------------------------
+
+#- shortcut to cfdem path
+alias cfdem='cd $CFDEM_PROJECT_DIR'
+
+#- shortcut to src path
+alias cfdemSrc='cd $CFDEM_SRC_DIR'
+
+#- shortcut to tutorial path
+alias cfdemTut='cd $CFDEM_TUT_DIR'
+
+#- shortcut to solver path
+alias cfdemSol='cd $CFDEM_SOLVER_DIR'
+
+#- shortcut to utilities path
+alias cfdemUt='cd $CFDEM_UT_DIR'
+
+#- shortcut to run path
+alias cfdemRun='cd $CFDEM_PROJECT_USER_DIR/run'
+
+#- shortcut to user solver path
+alias cfdemUsrSol='cd $CFDEM_PROJECT_USER_DIR/applications/solvers'
+
+#- shortcut to documentation path
+alias cfdemDoc='cd $CFDEM_DOC_DIR'
+
+#- shortcut to open the doxygen with firefox
+alias cfdemDox='firefox $CFDEM_DOC_DIR/doxygen/html/index.html'
+
+#- shortcut to LIGGGHTS path
+alias cfdemLIG='cd $CFDEM_LIGGGHTS_SRC_DIR'
+
+#- shortcut to system test
+alias cfdemSysTest='bash $CFDEM_PROJECT_DIR/etc/cfdemSystemTest.sh'
+
+#- shortcut to pull LIGGGHTS
+alias cfdemPullLIG='bash $CFDEM_PROJECT_DIR/etc/pullLIGGGHTS.sh'
+
+#- shortcut to pull CFDEMcoupling
+alias cfdemPullCFDEMcoupling='bash $CFDEM_PROJECT_DIR/etc/pullCFDEMcoupling.sh'
+
+#- shortcut to clean CFDEM
+alias cfdemCleanCFDEM='bash $CFDEM_PROJECT_DIR/etc/cleanCFDEMcoupling.sh'
+
+#- shortcut to compile LIGGGHTS + sublibraries
+alias cfdemCompLIG='bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS.sh'
+
+#- shortcut to compile CFDEMcoupling +LIGGGHTS
+alias cfdemCompCFDEMall='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_all.sh'
+
+#- shortcut to compile CFDEMcoupling (src+solvers)
+alias cfdemCompCFDEM='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling.sh'
+
+#- shortcut to compile CFDEMcoupling src
+alias cfdemCompCFDEMsrc='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_src.sh'
+
+#- shortcut to compile CFDEMcoupling solvers
+alias cfdemCompCFDEMsol='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_sol.sh'
+
+#- shortcut to compile CFDEMcoupling utilities
+alias cfdemCompCFDEMuti='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh'
+
+#- shortcut to test basic tutorials
+alias cfdemTestTUT='bash $CFDEM_PROJECT_DIR/etc/testTutorials.sh'
+
+#- shortcut to visualize the clock model data
+alias vizClock='python $CFDEM_UT_DIR/vizClock/matPlot.py'
+
+#- recursive touch of current directory
+alias touchRec='find ./* -exec touch {} \;'
+
+#- shortcut to run liggghts in serial
+cfdemLiggghts() { $CFDEM_LIGGGHTS_BIN_DIR/liggghts < $1; }
+export -f cfdemLiggghts
+
+#- shortcut to run liggghts in parallel
+cfdemLiggghtsPar() { mpirun -np $2 $CFDEM_LIGGGHTS_BIN_DIR/liggghts < $1; }
+export -f cfdemLiggghtsPar
+
+#- shortcut to open files including a pattern
+cfdemGrep() { grep -rl "$1" ./* | xargs gedit; }
+export -f cfdemGrep
+
+#- shortcut lo list files in a directory
+#cfdemListFiles() { find $1 | sed s:""$1"":: > listOfFiles.txt; }   #leave out the dir iteslf in list
+cfdemListFiles() { find $1 > listOfFiles.txt; }                     #keep the dir in list
+export -f cfdemListFiles
+
+#- check if the user directory exists
+if [ -d "$CFDEM_PROJECT_USER_DIR" ]; then
+    :
+else
+        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
+        read YN
+        if [ $YN == "y" ]; then
+            mkdir -p $CFDEM_PROJECT_USER_DIR
+            cd $CFDEM_PROJECT_USER_DIR
+            mkdir run
+            mkdir -p log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+            mkdir -p applications/solvers
+        else
+            echo "aborted by user."
+            #exit
+        fi
+fi

etc/cfdemSystemTest.sh

@@ -1,20 +1,25 @@
 #!/bin/bash
 
 #===================================================================#
-# sytsem settings test routine for cfdem project 
+# sytsem settings test routine for cfdem project
 # Christoph Goniva - May. 2011, DCS Computing GmbH
 #===================================================================#
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 #- show gcc settings
 checkGPP="true"
 
+#- sys check for add on
+checkAddOn="false"
+
 #- system settings
-echo "*******************"
-echo "system settings:"
-echo "*******************"
+printHeader
+
+echo "*********************************"
+echo "CFDEM(R)coupling system settings:"
+echo "*********************************"
 
 echo "CFDEM_VERSION=$CFDEM_VERSION"
 echo "couple to OF_VERSION=$WM_PROJECT_VERSION"
@@ -28,16 +33,17 @@ checkDirComment "$CFDEM_SRC_DIR" '$CFDEM_SRC_DIR' "yes"
 checkDirComment "$CFDEM_SOLVER_DIR" '$CFDEM_SOLVER_DIR' "yes"
 checkDirComment "$CFDEM_TUT_DIR" '$CFDEM_TUT_DIR' "yes"
 checkDirComment "$CFDEM_LIGGGHTS_SRC_DIR" '$CFDEM_LIGGGHTS_SRC_DIR' "yes"
+checkEnvComment "$CFDEM_LIGGGHTS_BIN_DIR" '$CFDEM_LIGGGHTS_BIN_DIR' "yes"
 checkDirComment "$CFDEM_LPP_DIR" '$CFDEM_LPP_DIR' "yes"
-checkDirComment "$CFDEM_PIZZA_DIR" '$CFDEM_PIZZA_DIR' "no"
+checkDirComment "$CFDEM_ADD_LIBS_DIR" '$CFDEM_ADD_LIBS_DIR' "yes"
 checkDirComment "$CFDEM_TEST_HARNESS_PATH" '$CFDEM_TEST_HARNESS_PATH' "no"
 echo ""
 
 echo "library names"
-echo '$CFDEM_LIGGGHTS_LIB_NAME = '"$CFDEM_LIGGGHTS_LIB_NAME"
 echo '$CFDEM_LIB_NAME = '"$CFDEM_LIB_NAME"
-echo '$LD_LIBRARY_PATH  = '"$LD_LIBRARY_PATH"
-echo '$WM_NCOMPPROCS  = '"$WM_NCOMPPROCS"
+echo '$PATH = '"$PATH"
+echo '$LD_LIBRARY_PATH = '"$LD_LIBRARY_PATH"
+echo '$WM_NCOMPPROCS = '"$WM_NCOMPPROCS"
 
 echo "*******************"
 
@@ -61,3 +67,16 @@ if [ $checkGPP == "true" ]
     mpirun --version
 fi
 
+if [ $checkAddOn == "true" ]
+  then
+    echo "**********************"
+    echo "additional packages..."
+
+    packageName=c3po
+    filePath=$CFDEM_SRC_DIR/$packageName
+    if [ $(checkDir $filePath) == "true" ]; then
+        source $filePath/etc/$packageName"SystemTest.sh"
+    else
+        echo "$packageName does not exist."
+    fi
+fi

etc/cleanCFDEMcoupling.sh

@@ -6,7 +6,7 @@
 #===================================================================#
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 #================================================================================#
 # clean src remove object files

etc/compileCFDEMcoupling.sh

@@ -6,26 +6,26 @@
 #===================================================================#
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #================================================================================#
 # compile src
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_src.sh
+bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_src.sh
 
 #================================================================================#
 # compile solvers
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_sol.sh
+bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_sol.sh
 
 #================================================================================#
 # compile utilities
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh
+bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh

etc/compileCFDEMcoupling_all.sh

@@ -7,26 +7,26 @@
 #===================================================================#
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #================================================================================#
 # compile LIGGGHTS src
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS.sh
+bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS.sh
 
 #================================================================================#
 # compile LIGGGHTS libraries
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh
+bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS_lib.sh
 
 #================================================================================#
 # compile CFDEMcoupling
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling.sh
+bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling.sh

etc/compileCFDEMcoupling_sol.sh

@@ -8,9 +8,9 @@
 whitelist="solver-list.txt"
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 logDir="log"
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #- remove old success/fail logs

etc/compileCFDEMcoupling_src.sh

@@ -7,20 +7,20 @@
 #===================================================================#
  
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 
 #================================================================================#
 # compile src
 #================================================================================#
-whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/library-list.txt"
+whitelist="$CFDEM_PROJECT_DIR/etc/library-list.txt"
 echo ""
 echo "Please provide the libraries to be compiled in the $CWD/$whitelist file."
 

etc/compileCFDEMcoupling_uti.sh

@@ -8,9 +8,9 @@
 whitelist="utilities-list.txt"
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 logDir="log"
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 CWD="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"

etc/compileLIGGGHTS.sh

@@ -6,12 +6,12 @@
 #===================================================================
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #================================================================================#
@@ -38,4 +38,4 @@ compileLIGGGHTS $logpath $logfileName $headerText
 #================================================================================#
 # compile LIGGGHTS libraries
 #================================================================================#
-bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh
+bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS_lib.sh

etc/compileLIGGGHTS_lib.sh

@@ -6,19 +6,19 @@
 #===================================================================#
  
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #================================================================================#
 # compile src
 #================================================================================#
-    whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/library-liggghts-list.txt"
+    whitelist="$CFDEM_PROJECT_DIR/etc/library-liggghts-list.txt"
     echo ""
     echo "Compiling sub-libraries of LIGGGHTS now..."
     echo "Please provide the libraries to be compiled in the $CWD/$whitelist file."

etc/cshrc

@@ -0,0 +1,289 @@
+#----------------------------------*-sh-*--------------------------------------
+# CFDEMcoupling
+#
+# Copyright 2009-2012 JKU Linz
+# Copyright 2012-2015 DCS Computing GmbH, Linz
+# Copyright 2015-     JKU Linz
+#------------------------------------------------------------------------------
+#
+# File
+#     etc/cshrc
+#
+# Description
+#     Startup file for CFDEMcoupling
+#     Sourced from ~/.cshrc
+#------------------------------------------------------------------------------
+
+setenv CFDEM_PROJECT CFDEM
+setenv CFDEM_VERSION 17.08
+
+################################################################################
+# USER EDITABLE PART: Changes made here may be lost with the next upgrade
+#
+# Please set to the appropriate path if the default is not correct.
+#
+# activate compatible OpenFOAM version
+. $HOME/OpenFOAM/OpenFOAM-4.x/etc/cshrc
+#
+# CFDEMcoupling
+setenv CFDEM_INST_DIR $HOME/$CFDEM_PROJECT
+setenv CFDEM_PROJECT_DIR $CFDEM_INST_DIR/CFDEMcoupling
+setenv CFDEM_PROJECT_USER_DIR $CFDEM_INST_DIR/$USER-$WM_PROJECT_VERSION
+#
+# LIGGGHTS
+setenv CFDEM_LIGGGHTS_INST_DIR $CFDEM_INST_DIR/LIGGGHTS
+#
+# LPP installation path
+setenv CFDEM_LPP_INST_DIR $CFDEM_INST_DIR/LPP
+#
+# Path to additional libraries
+setenv CFDEM_ADD_LIBS_DIR
+#
+# END OF (NORMAL) USER EDITABLE PART
+################################################################################
+
+#- The old dirs to be cleaned from the environment variables
+set cfdemOldDirs=
+if ( $?CFDEM_LIGGGHTS_BIN_DIR ) then
+   set cfdemOldDirs="$cfdemOldDirs $CFDEM_LIGGGHTS_BIN_DIR"
+endif
+if ( $?CFDEM_APP_DIR ) then
+   set cfdemOldDirs="$cfdemOldDirs $CFDEM_APP_DIR"
+endif
+if ( $?CFDEM_LIB_DIR ) then
+   set cfdemOldDirs="$cfdemOldDirs $CFDEM_LIB_DIR"
+endif
+
+setenv CFDEM_LIGGGHTS_SRC_DIR $CFDEM_LIGGGHTS_INST_DIR/src
+setenv CFDEM_LIGGGHTS_BIN_DIR $CFDEM_LIGGGHTS_INST_DIR/src-build
+
+#------------------------------------------------------------------------------
+# Source files, possibly with some verbosity
+#------------------------------------------------------------------------------
+alias _foamSource 'if ($?FOAM_VERBOSE && $?prompt) echo "Sourcing: \!*"; if (\!* != "") source \!*'
+
+#------------------------------------------------------------------------------
+# Evaluate command-line parameters
+# these can be used to set/unset values
+#------------------------------------------------------------------------------
+while ( $#argv > 0 )
+    switch ($argv[1])
+    case -*:
+        # stray option (not meant for us here) -> get out
+        break
+        breaksw
+    case *=:
+        # name=       -> unsetenv name
+        if ($?FOAM_VERBOSE && $?prompt) echo "unsetenv $argv[1]:s/=//"
+        eval "unsetenv $argv[1]:s/=//"
+        breaksw
+    case *=*:
+        # name=value  -> setenv name value
+        if ($?FOAM_VERBOSE && $?prompt) echo "setenv $argv[1]:s/=/ /"
+        eval "setenv $argv[1]:s/=/ /"
+        breaksw
+    default:
+        # filename: source it
+        if ( -f "$1" ) then
+            _foamSource "$1"
+        else
+            _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile -silent "$1"`
+        endif
+        breaksw
+    endsw
+    shift
+end
+
+#------------------------------------------------------------------------------
+
+setenv CFDEM_SRC_DIR $CFDEM_PROJECT_DIR/src
+setenv CFDEM_SOLVER_DIR $CFDEM_PROJECT_DIR/applications/solvers
+setenv CFDEM_UT_DIR $CFDEM_PROJECT_DIR/applications/utilities
+setenv CFDEM_DOC_DIR $CFDEM_PROJECT_DIR/doc
+setenv CFDEM_TUT_DIR $CFDEM_PROJECT_DIR/tutorials
+setenv CFDEM_LPP_DIR $CFDEM_LPP_INST_DIR/src
+
+#------------------------------------------------------------------------------
+
+#- CFDEM lib name
+setenv CFDEM_LIB_NAME lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#- CFDEM compressible lib name
+setenv CFDEM_LIB_COMP_NAME lagrangianCFDEMcomp-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#- Check if additional libraries should be compiled together with solvers
+if ( ! ($?CFDEM_ADD_LIBS_DIR) ) then
+    setenv CFDEM_ADD_LIBS_DIR $CFDEM_PROJECT_DIR/etc
+else
+    echo "using CFDEM_ADD_LIBS_DIR=$CFDEM_ADD_LIBS_DIR defined by user."
+endif
+
+#------------------------------------------------------------------------------
+
+#- LMP Many2Many lib path and makefile
+setenv CFDEM_Many2ManyLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
+setenv CFDEM_Many2ManyLIB_MAKEFILENAME fedora_fpic
+
+#------------------------------------------------------------------------------
+
+#- path to test harness
+setenv CFDEM_TEST_HARNESS_PATH $CFDEM_PROJECT_USER_DIR/log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+
+#- path to libraries
+setenv CFDEM_LIB_DIR $CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/lib
+
+#- path to apps
+setenv CFDEM_APP_DIR $CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/bin
+
+#- create directories
+mkdir -p $CFDEM_LIB_DIR
+mkdir -p $CFDEM_APP_DIR
+
+#- path to OF version flag file
+setenv CFDEM_OFVERSION_DIR $CFDEM_PROJECT_DIR/etc/OFversion
+
+#------------------------------------------------------------------------------
+# Clean standard environment variables (PATH, LD_LIBRARY_PATH)
+#------------------------------------------------------------------------------
+set foamClean=$WM_PROJECT_DIR/bin/foamCleanPath
+
+#- prevent local variables from shadowing setenv variables
+unset PATH LD_LIBRARY_PATH
+
+if (! $?LD_LIBRARY_PATH ) setenv LD_LIBRARY_PATH ''
+
+#- Clean PATH
+set cleaned=`$foamClean "$PATH" "$cfdemOldDirs"`
+if ( $status == 0 ) setenv PATH $cleaned
+
+#- Clean LD_LIBRARY_PATH
+set cleaned=`$foamClean "$LD_LIBRARY_PATH" "$cfdemOldDirs"`
+if ( $status == 0 ) setenv LD_LIBRARY_PATH $cleaned
+
+#- add binary directories to $PATH
+_foamAddPath $CFDEM_APP_DIR:$CFDEM_LIGGGHTS_BIN_DIR
+_foamAddLib $CFDEM_LIB_DIR
+
+#------------------------------------------------------------------------------
+# Cleanup environment:
+#------------------------------------------------------------------------------
+unset cleaned foamClean cfdemOldDirs
+unalias _foamSource
+
+#------------------------------------------------------------------------------
+#- settings for lpp postproc tool
+#------------------------------------------------------------------------------
+
+#- nr of procs for lpp tool
+setenv CFDEM_LPP_NPROCS 4
+
+#- nr of procs for lpp tool
+setenv CFDEM_LPP_CHUNKSIZE 1
+
+#- shortcut to run lpp
+alias lpp 'python -i $CFDEM_LPP_DIR/lpp.py --cpunum $CFDEM_LPP_NPROCS --chunksize $CFDEM_LPP_CHUNKSIZE \!:1'
+
+#------------------------------------------------------------------------------
+# aliases for easy navigation (no changes necessary)
+#------------------------------------------------------------------------------
+
+#- shortcut to cfdem path
+alias cfdem 'cd $CFDEM_PROJECT_DIR'
+
+#- shortcut to src path
+alias cfdemSrc 'cd $CFDEM_SRC_DIR'
+
+#- shortcut to tutorial path
+alias cfdemTut 'cd $CFDEM_TUT_DIR'
+
+#- shortcut to solver path
+alias cfdemSol 'cd $CFDEM_SOLVER_DIR'
+
+#- shortcut to utilities path
+alias cfdemUt 'cd $CFDEM_UT_DIR'
+
+#- shortcut to run path
+alias cfdemRun 'cd $CFDEM_PROJECT_USER_DIR/run'
+
+#- shortcut to user solver path
+alias cfdemUsrSol 'cd $CFDEM_PROJECT_USER_DIR/applications/solvers'
+
+#- shortcut to documentation path
+alias cfdemDoc 'cd $CFDEM_DOC_DIR'
+
+#- shortcut to open the doxygen with firefox
+alias cfdemDox 'firefox $CFDEM_DOC_DIR/doxygen/html/index.html'
+
+#- shortcut to LIGGGHTS path
+alias cfdemLIG 'cd $CFDEM_LIGGGHTS_SRC_DIR'
+
+#- shortcut to system test
+alias cfdemSysTest 'bash $CFDEM_PROJECT_DIR/etc/cfdemSystemTest.sh'
+
+#- shortcut to pull LIGGGHTS
+alias cfdemPullLIG 'bash $CFDEM_PROJECT_DIR/etc/pullLIGGGHTS.sh'
+
+#- shortcut to pull CFDEMcoupling
+alias cfdemPullCFDEMcoupling 'bash $CFDEM_PROJECT_DIR/etc/pullCFDEMcoupling.sh'
+
+#- shortcut to clean CFDEM
+alias cfdemCleanCFDEM 'bash $CFDEM_PROJECT_DIR/etc/cleanCFDEMcoupling.sh'
+
+#- shortcut to compile LIGGGHTS + sublibraries
+alias cfdemCompLIG 'bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS.sh'
+
+#- shortcut to compile CFDEMcoupling +LIGGGHTS
+alias cfdemCompCFDEMall 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_all.sh'
+
+#- shortcut to compile CFDEMcoupling (src+solvers)
+alias cfdemCompCFDEM 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling.sh'
+
+#- shortcut to compile CFDEMcoupling src
+alias cfdemCompCFDEMsrc 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_src.sh'
+
+#- shortcut to compile CFDEMcoupling solvers
+alias cfdemCompCFDEMsol 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_sol.sh'
+
+#- shortcut to compile CFDEMcoupling utilities
+alias cfdemCompCFDEMuti 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh'
+
+#- shortcut to test basic tutorials
+alias cfdemTestTUT 'bash $CFDEM_PROJECT_DIR/etc/testTutorials.sh'
+
+#- shortcut to visualize the clock model data
+alias vizClock 'python $CFDEM_UT_DIR/vizClock/matPlot.py'
+
+#- shortcut to run liggghts in serial
+alias cfdemLiggghts '$CFDEM_LIGGGHTS_BIN_DIR/liggghts'
+
+#- shortcut to run liggghts in parallel (no fcts in csh...)
+alias cfdemLiggghtsPar 'echo "mpirun -np xx -machinefile mynodes $CFDEM_LIGGGHTS_BIN_DIR/liggghts < in.liggghts_init"'
+
+#- shortcut to run liggghts in parallel
+# unfortunately no functions available in csh
+
+#- shortcut to run lpp
+alias lpp 'python -i $CFDEM_LPP_DIR/lpp.py \!:1'
+
+#- check if the user directory exists
+if ( -d "$CFDEM_PROJECT_USER_DIR" ) then
+    :
+else
+    if ( -d "$CFDEM_PROJECT_DIR" ) then
+        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
+        set YN=$<
+        if ( $YN == "y" ) then
+            mkdir -p $CFDEM_PROJECT_USER_DIR
+            cd $CFDEM_PROJECT_USER_DIR
+            mkdir run
+            mkdir -p log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
+            mkdir -p applications/solvers
+        else
+            echo "aborted by user."
+            exit
+        endif
+    else
+        echo "error in CFDEMcoupling's cshrc."
+        exit
+    endif
+endif

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 "Please make sure to have the compressible libraries first in the library-list.txt!"
+                echo "Compiling a compressible library - so doing a wrmdep -a of incompressible library first."
+                echo "Please make sure to have the incompressible 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."
+                echo "Compiling an 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
     
@@ -161,7 +161,8 @@ compileLIGGGHTS()
     rm $logpath/$logfileName
 
     #- change path
-    cd $CFDEM_LIGGGHTS_SRC_DIR
+    mkdir -p $CFDEM_LIGGGHTS_BIN_DIR
+    cd $CFDEM_LIGGGHTS_BIN_DIR
 
     #- header
     echo 2>&1 | tee -a $logpath/$logfileName
@@ -176,20 +177,18 @@ compileLIGGGHTS()
     if [[ $clean == "false" ]]; then
         echo "not cleaning LIGGGHTS"
     else
-        rm $CFDEM_LIGGGHTS_SRC_DIR/"lmp_"$CFDEM_LIGGGHTS_MAKEFILE_NAME
-        rm $CFDEM_LIGGGHTS_SRC_DIR/"lib"$CFDEM_LIGGGHTS_LIB_NAME".a"
-        make clean-all 2>&1 | tee -a $logpath/$logfileName
+        make clean 2>&1 | tee -a $logpath/$logfileName
+        rm CMakeCache.txt 2>&1 | tee -a $logpath/$logfileName
         echo "cleaning LIGGGHTS"
     fi
+    cmake $CFDEM_LIGGGHTS_SRC_DIR
     if [[ $WM_NCOMPPROCS == "" ]]; then
         echo "compiling LIGGGHTS on one CPU"
-        make $CFDEM_LIGGGHTS_MAKEFILE_NAME 2>&1 | tee -a $logpath/$logfileName
+        make 2>&1 | tee -a $logpath/$logfileName
     else
         echo "compiling LIGGGHTS on $WM_NCOMPPROCS CPUs"
-        make $CFDEM_LIGGGHTS_MAKEFILE_NAME -j $WM_NCOMPPROCS  2>&1 | tee -a $logpath/$logfileName
+        make -j $WM_NCOMPPROCS  2>&1 | tee -a $logpath/$logfileName
     fi
-    make makelib 2>&1 | tee -a $logpath/$logfileName
-    make -f Makefile.lib $CFDEM_LIGGGHTS_MAKEFILE_NAME 2>&1 | tee -a $logpath/$logfileName
 }
 
 #==================================#
@@ -250,7 +249,7 @@ cleanCFDEM()
 
     #**********************************************
     #cleaning libraries
-    whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/library-list.txt"
+    whitelist="$CFDEM_PROJECT_DIR/etc/library-list.txt"
     echo ""
     echo "Please provide the libraries to be cleaned in the $CWD/$whitelist file."
 
@@ -285,7 +284,7 @@ cleanCFDEM()
 
             cd  $path
             echo "cleaning library $PWD"
-            rmdepall
+            wrmdep -a
             wclean    
             rm -r ./Make/linux*
             rm -r ./lnInclude
@@ -304,7 +303,7 @@ cleanCFDEM()
 
     #**********************************************
     #cleaning solvers
-    whitelist="$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/solver-list.txt"
+    whitelist="$CFDEM_PROJECT_DIR/etc/solver-list.txt"
     echo ""
     echo "Please provide the solvers to be cleaned in the $CWD/$whitelist file."
 
@@ -339,7 +338,7 @@ cleanCFDEM()
 
             cd  $path            
             echo "cleaning solver $PWD"
-            rmdepall
+            wrmdep -a
             wclean    
     done
 }
@@ -405,7 +404,7 @@ DEMrun()
     echo 2>&1 | tee -a $logpath/$logfileName
 
     #- run applictaion
-    $debugMode $CFDEM_LIGGGHTS_SRC_DIR/$CFDEM_LIGGGHTS_LIB_NAME < $solverName 2>&1 | tee -a $logpath/$logfileName
+    $debugMode $CFDEM_LIGGGHTS_BIN_DIR/liggghts -in $solverName 2>&1 | tee -a $logpath/$logfileName
 
     #- keep terminal open (if started in new terminal)
     #read
@@ -455,9 +454,9 @@ parDEMrun()
 
     #- run applictaion
     if [ $machineFileName == "none" ]; then
-        mpirun -np $nrProcs $debugMode $CFDEM_LIGGGHTS_SRC_DIR/$CFDEM_LIGGGHTS_LIB_NAME < $solverName 2>&1 | tee -a $logpath/$logfileName
+        mpirun -np $nrProcs $debugMode $CFDEM_LIGGGHTS_BIN_DIR/liggghts -in $solverName 2>&1 | tee -a $logpath/$logfileName
     else
-        mpirun -machinefile $machineFileName -np $nrProcs $debugMode $CFDEM_LIGGGHTS_SRC_DIR/$CFDEM_LIGGGHTS_LIB_NAME < $solverName 2>&1 | tee -a $logpath/$logfileName
+        mpirun -machinefile $machineFileName -np $nrProcs $debugMode $CFDEM_LIGGGHTS_BIN_DIR/liggghts -in $solverName 2>&1 | tee -a $logpath/$logfileName
     fi
 
     #- keep terminal open (if started in new terminal)
@@ -931,6 +930,52 @@ checkDirComment()
     fi
 }
 
+#========================================#
+#- function to check if a variable exits
+checkEnv()
+{
+    #--------------------------------------------------------------------------------#
+    #- define variables
+    var="$1"
+    #--------------------------------------------------------------------------------#
+    if [[ $var == "" ]]; then
+        echo "false"
+    else
+        echo "true"
+    fi
+}
+
+#========================================#
+#- function to check if a variable exits
+checkEnvComment()
+{
+    #--------------------------------------------------------------------------------#
+    #- define variables
+    var="$1"
+    varName="$2"
+    critical="$3"
+    #--------------------------------------------------------------------------------#
+    if [ $(checkEnv $var) == "true" ]; then
+         echo "valid:yes critical:$critical - $varName = $var" 
+    else
+        echo "valid:NO  critical:$critical - $varName = $var variable not set!" 
+    fi
+}
+
+#========================================#
+#- function to print a header to terminal
+printHeader()
+{
+    echo ""
+    echo "*********************************************"
+    echo "* C F D E M (R) c o u p l i n g             *"
+    echo "*                                           *"
+    echo "* by DCS Computing GmbH                     *"
+    echo "* www.dcs-computing.com                     *"
+    echo "*********************************************"
+    echo "" 
+}
+
 #========================================#
 #- track memory usage
 trackMem()

etc/library-list.txt

@@ -1,2 +1,3 @@
 lagrangian/cfdemParticle/dir
+lagrangian/cfdemParticleComp/dir
 finiteVolume/dir

etc/pullLIGGGHTS.sh

@@ -6,12 +6,12 @@
 #===================================================================
 
 #- include functions
-source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/functions.sh
+source $CFDEM_PROJECT_DIR/etc/functions.sh
 
 NOW="$(date +"%Y-%m-%d-%H:%M")"
 logDir="log"
 
-cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
+cd $CFDEM_PROJECT_DIR/etc
 mkdir -p $logDir
 
 #--------------------------------------------------------------------------------#

etc/solver-list.txt

@@ -1,4 +1,5 @@
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
+cfdemSolverRhoPimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir

src/finiteVolume/Make/options

@@ -2,6 +2,9 @@ EXE_INC = \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -Wno-old-style-cast \
+    -Wno-unused-result \
+    -Wno-literal-suffix
 
 LIB_LIBS = \
     -lOpenFOAM \

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"),
@@ -195,6 +195,7 @@ template<class Type>
 void uniformFixedValueTubeFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
+    this->writeEntry("value", os);
     uniformValue_->writeData(os);
     os.writeKeyword("tubeLength") << tubeLength_ << token::END_STATEMENT << nl;
     os.writeKeyword("tubeDiameter") << tubeDiameter_ << token::END_STATEMENT << nl;

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

@@ -41,7 +41,6 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 :
     fixedValueFvPatchField<Type>(p, iF),
     uniformValue_(),
-    //voidfractionFieldName_(dict.lookupOrDefault<word>("voidfraction", "voidfraction"))
     voidfractionFieldName_("voidfraction")
 {}
 
@@ -56,7 +55,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(ptf.uniformValue_().clone().ptr()),
+    uniformValue_(ptf.uniformValue_, false),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();
@@ -73,7 +72,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();
@@ -88,7 +87,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(ptf),
-    uniformValue_(ptf.uniformValue_().clone().ptr()),
+    uniformValue_(ptf.uniformValue_, false),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();
@@ -104,7 +103,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(ptf, iF),
-    uniformValue_(ptf.uniformValue_().clone().ptr()),
+    uniformValue_(ptf.uniformValue_, false),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();
@@ -146,6 +145,7 @@ template<class Type>
 void uniformFixedValueVoidfractionFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
+    this->writeEntry("value", os);
     uniformValue_->writeData(os);
 }
 

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/.gitignore

@@ -0,0 +1 @@
+lnInclude

src/lagrangian/cfdemParticle/Make/files

@@ -1,8 +1,11 @@
 cfdemCloud = cfdemCloud
 cfdTools = cfdTools
+energyModels = subModels/energyModel
 forceModels = subModels/forceModel
 forceSubModels = subModels/forceModel/forceSubModels
 forceModelsMS = subModels/forceModelMS
+thermCondModels = subModels/thermCondModel
+chemistryModels = subModels/chemistryModel
 IOModels = subModels/IOModel
 voidFractionModels = subModels/voidFractionModel
 locateModels = subModels/locateModel
@@ -15,14 +18,28 @@ clockModels = subModels/clockModel
 liggghtsCommandModels = subModels/liggghtsCommandModel
 smoothingModels = subModels/smoothingModel
 probeModels = subModels/probeModel
+otherForceModels = subModels/otherForceModel
 
 $(cfdemCloud)/cfdemCloud.C
+$(cfdemCloud)/cfdemCloudIO.C
 derived/cfdemCloudIB/cfdemCloudIB.C
 derived/cfdemCloudMS/cfdemCloudMS.C
+derived/cfdemCloudEnergy/cfdemCloudEnergy.C
 
 $(cfdTools)/global.C
 $(cfdTools)/newGlobal.C
 
+$(energyModels)/energyModel/energyModel.C
+$(energyModels)/energyModel/newEnergyModel.C
+$(energyModels)/heatTransferGunn/heatTransferGunn.C
+$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C
+$(energyModels)/reactionHeat/reactionHeat.C
+
+$(thermCondModels)/thermCondModel/thermCondModel.C
+$(thermCondModels)/thermCondModel/newThermCondModel.C
+$(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
+$(thermCondModels)/noTherm/noThermCond.C
+
 $(forceModels)/forceModel/forceModel.C
 $(forceModels)/forceModel/newForceModel.C
 $(forceModels)/noDrag/noDrag.C
@@ -45,6 +62,13 @@ $(forceModels)/MeiLift/MeiLift.C
 $(forceModels)/particleCellVolume/particleCellVolume.C
 $(forceModels)/fieldTimeAverage/fieldTimeAverage.C
 $(forceModels)/volWeightedAverage/volWeightedAverage.C
+$(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/dSauter/dSauter.C
+$(forceModels)/Fines/Fines.C
+$(forceModels)/Fines/FinesFields.C
+$(forceModels)/Fines/FanningDynFines.C
+$(forceModels)/Fines/ErgunStatFines.C
+$(forceModels)/granKineticEnergy/granKineticEnergy.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
@@ -53,6 +77,19 @@ $(forceModelsMS)/DiFeliceDragMS/DiFeliceDragMS.C
 $(forceSubModels)/forceSubModel/newForceSubModel.C
 $(forceSubModels)/forceSubModel/forceSubModel.C
 $(forceSubModels)/ImEx/ImEx.C
+$(forceSubModels)/ScaleForce/ScaleForce.C
+$(forceSubModels)/scaleForceBoundary/scaleForceBoundary.C
+
+$(otherForceModels)/otherForceModel/otherForceModel.C
+$(otherForceModels)/otherForceModel/newOtherForceModel.C
+$(otherForceModels)/gravity/gravity.C
+$(otherForceModels)/weightSecondaryPhase/weightSecondaryPhase.C
+$(otherForceModels)/expParticleForces/expParticleForces.C
+
+$(chemistryModels)/chemistryModel/chemistryModel.C
+$(chemistryModels)/chemistryModel/newChemistryModel.C
+$(chemistryModels)/species/species.C
+$(chemistryModels)/noChemistry/noChemistry.C
 
 $(probeModels)/probeModel/probeModel.C
 $(probeModels)/probeModel/newProbeModel.C
@@ -74,6 +111,7 @@ $(voidFractionModels)/dividedVoidFractionMS/dividedVoidFractionMS.C
 $(voidFractionModels)/bigParticleVoidFraction/bigParticleVoidFraction.C
 $(voidFractionModels)/GaussVoidFraction/GaussVoidFraction.C
 $(voidFractionModels)/IBVoidFraction/IBVoidFraction.C
+$(voidFractionModels)/trilinearVoidFraction/trilinearVoidFraction.C
 
 $(locateModels)/locateModel/locateModel.C
 $(locateModels)/locateModel/newLocateModel.C
@@ -83,7 +121,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
@@ -104,7 +141,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

@@ -1,34 +1,39 @@
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
 sinclude $(RULES)/mplib$(WM_MPLIB)
 
-GIT_VERSION := $(shell git describe --abbrev=4 --dirty --always --tags)
+GIT_VERSION := $(shell git describe --dirty --always --tags)
 PFLAGS+= -DGITVERSION=\"$(GIT_VERSION)\"
 
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
      $(PFLAGS) \
      $(PINC) \
-    -I ./cfdemParticle \
+    -I$(CFDEM_OFVERSION_DIR) \
     -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 \
     -I$(CFDEM_LIGGGHTS_SRC_DIR) \
     -I$(CFDEM_M2MLIB_PATH) \
     -I$(CFDEM_SRC_DIR)/cfdTools \
+    -Wno-old-style-cast \
+    -Wno-unused-result \
+    -Wno-literal-suffix
+
 
 LIB_LIBS = \
      $(PLIBS) \
     -L$(CFDEM_LIB_DIR) \
     -lfiniteVolume \
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lincompressibleTurbulenceModels \
     -lmeshTools \
     -llagrangian \
     -lmpi_cxx \
-    -L$(CFDEM_LIGGGHTS_SRC_DIR) \
-    -Wl,--whole-archive -l$(CFDEM_LIGGGHTS_LIB_NAME) -Wl,--no-whole-archive \
+    -Wl,-rpath,$(CFDEM_LIGGGHTS_BIN_DIR) \
+    -L$(CFDEM_LIGGGHTS_BIN_DIR) \
+    -lliggghts \
     -L$(CFDEM_Many2ManyLIB_PATH) \
     -lcoupleMany2Many

src/lagrangian/cfdemParticle/cfdTools/compressibleContinuityErrsPU.H

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-     DCS Computing GmbH,Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    continuityErrs
+
+Description
+    Calculates and prints the continuity errors.
+    The code is an evolution of compressibleContinuityErrs.H in OpenFOAM(R) 2.3.x,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+{
+    dimensionedScalar totalMass = fvc::domainIntegrate(rho*voidfraction);
+
+    scalar sumLocalContErr =
+        (fvc::domainIntegrate(mag(rho - thermo.rho())*voidfraction)/totalMass).value();
+
+    scalar globalContErr =
+        (fvc::domainIntegrate((rho - thermo.rho())*voidfraction)/totalMass).value();
+
+    cumulativeContErr += globalContErr;
+
+    Info<< "time step continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}
+
+// ************************************************************************* //

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/cfdTools/debugInfo.H

@@ -1,7 +1,7 @@
-
+{
   //========================================================================//
-  scalar countCell=0;              // number of cells touched by particles
-  int points=0;                    // number of particles and sub-points
+  label countCell=0;               // number of cells touched by particles
+  label points=0;                  // number of particles and sub-points
   scalar totalParticleWeights=0;   // total weight of all particles and sub-points
   vector totalForce_array(0,0,0);  // total force on particles based on particle array
   vector totalForce_field(0,0,0);  // forceField of forceM(), used to calc Ksl
@@ -23,8 +23,8 @@
           meanUs_array[i] += velocities_[index][i];
       }
       meanR_array += radii_[index][0];
-      particleVolume_radius += radii_[index][0]*radii_[index][0]*radii_[index][0]*4./3.*3.1415;
-      
+      particleVolume_radius += radii_[index][0]*radii_[index][0]*radii_[index][0]*(4./3.)*M_PI;
+
       // loop subCells
       for(int subCell=0;subCell<voidFractionM().cellsPerParticle()[index][0];subCell++)
       {
@@ -63,17 +63,17 @@
   if(countCell>0)
   {
       meanAlpha_field /= countCell;
-      meanU_field /= countCell; 
+      meanU_field /= countCell;
       meanUs_field /= countCell;
   }
   else
   {
       meanAlpha_field = 0;
-      meanU_field = vector(0,0,0); 
+      meanU_field = vector(0,0,0);
       meanUs_field = vector(0,0,0);
   }
-  meanUs_array /= numberOfParticles()+SMALL;   
-  meanR_array /= numberOfParticles()+SMALL;   
+  meanUs_array /= numberOfParticles()+SMALL;
+  meanR_array /= numberOfParticles()+SMALL;
 
   Info <<"=============================================================================" << endl;
   Info << "Debug Info, only serial and not tested!" << endl;
@@ -93,3 +93,4 @@
   Info <<"meanR_array = "<< meanR_array << endl;
   Info <<"=============================================================================" << endl;
   Info << endl;
+}

src/lagrangian/cfdemParticle/cfdTools/global.C

@@ -47,8 +47,8 @@ defineRunTimeSelectionTable(global, dictionary);
 void global::info()
 {
     Info << "\nYou are currently using:" << endl;
-    Info << "OF version: " << FOAMversion << endl;
-    Info << "OF build: " << FOAMbuild << endl;
+    Info << "OpenFOAM version: " << FOAMversion << endl;
+    Info << "OpenFOAM build: " << FOAMbuild << endl;
     Info << "CFDEM build: " << CFDEMversion << "\n" << endl;
 }
 

src/lagrangian/cfdemParticle/cfdTools/mathExtra.H

@@ -27,7 +27,7 @@ License
 Description
     This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
     and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
-    
+
     Copyright of this contribution:
     Copyright 2014-     TU Graz, IPPT
 ------------------------------------------------------------------------- */
@@ -36,15 +36,37 @@ Description
 #define CFDEM_MATH_EXTRA_H
 
 #include <vector>
-//#include "math.h"
-#include "stdio.h"
-#include "string.h"
-#include "error.h"
-#include "ctype.h"
+#include <stdio.h>
+#include <string.h>
+#include <error.h>
+#include <ctype.h>
+#include "mathematicalConstants.H"
 
 #define TOLERANCE_ORTHO 1e-10
 
-namespace MathExtra 
+namespace Foam
+{
+namespace constant
+{
+namespace mathematical
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    const scalar piByFour(0.25*pi);
+    const scalar fourPiByThree(4.0*pi/3.0);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace mathematical
+} // End namespace constant
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+namespace MathExtra
 {
 
 //  inline void     outerProduct(double *vec1, double *vec2, double **m);
@@ -58,11 +80,11 @@ inline void outerProduct(double *vec1, double *vec2, double **m)
 {
   int i, j;
   //debug output
-  for( i = 0; i < 3; ++i )
+//  for( i = 0; i < 3; ++i )
 //    printf("OUTER PRODUCT: Input: vec1 element %d = %g", i, vec1[i]);
-  for( i = 0; i < 3; ++i )
+//  for( i = 0; i < 3; ++i )
 //    printf("OUTER PRODUCT: Input: vec2 element %d=%g", i, vec2[i]);
-  
+
   //calculation
   for( i = 0; i < 3; ++i )
     for( j = 0; j < 3; ++j )
@@ -70,41 +92,33 @@ inline void outerProduct(double *vec1, double *vec2, double **m)
       m[i][j] = vec1[i] * vec2[j];
       printf("OUTER PRODUCT: Result: m[%d][%d]=%g", i, j, m[i][j]);
     }
-    
-    
 }
 //--------------------------------------------------------------------
 // Compute the major, minor axis and eccentricity parameters of a prolate spheroid
-inline bool spheroidGeometry(double radius, double aspectRatio,  //inputs
-                                 double& ai,    double& bi,              //outputs
-                                 double& ei,    double& Le              //outputs
-                                ) //
+inline bool spheroidGeometry(double radius, double aspectRatio, //inputs
+                             double& ai,    double& bi,         //outputs
+                             double& ei,    double& Le          //outputs
+                            )
 {
     //INPUT
     //  radius      ...volume-equivalent radius of the spheroid
     //  aspectRatio ...major/minor aspect ratio
-    
+
     //OUTPUT
-    //  ai  ... 
-    //  bi  ... 
-    //  ei  ... 
-    //  Le  ... 
-  
+    //  ai  ...
+    //  bi  ...
+    //  ei  ...
+    //  Le  ...
+
     if(radius<=0.0)     //avoid troubles in case radius is 0 or negative
         return false;
 
- 	ai = radius * std::pow(aspectRatio*aspectRatio,0.33333333333333333333333);
+    ai = radius * std::pow(aspectRatio*aspectRatio,0.33333333333333333333333);
     bi = ai / aspectRatio;
- 	ei = std::sqrt( 
- 	                1.0
-                   - 1.0 / (aspectRatio*aspectRatio)
- 	              );
- 	Le = std::log(
- 		           (1.0+ei)
- 		          /(1.0-ei)
-                 );
+    ei = std::sqrt( 1.0 - 1.0 / (aspectRatio*aspectRatio) );
+    Le = std::log( (1.0+ei) / (1.0-ei) );
 
-	return true; 		    
+    return true;
 }
 
 //--------------------------------------------------------------------
@@ -117,17 +131,17 @@ inline double Pi()
 
 //--------------------------------------------------------------------
 // Compute the major, minor axis and eccentricity parameters of a prolate spheroid
-inline bool spheroidGeometry2(double radius, double aspectRatio,  //inputs
-                                 double& ai,     double& bi,           //outputs
-                                 double& XAe,    double& YAe,         //outputs
-                                 double& XCe,    double& YCe,         //outputs
+inline bool spheroidGeometry2(double radius,  double aspectRatio, //inputs
+                                 double& ai,  double& bi,         //outputs
+                                 double& XAe, double& YAe,        //outputs
+                                 double& XCe, double& YCe,        //outputs
                                  double& YHe
-                                ) //
+                                )
 {
     //INPUT
     //  radius      ...volume-equivalent radius of the spheroid
     //  aspectRatio ...major/minor aspect ratio
-    
+
     //OUTPUT
     //  XAe  ...Eccentricity dependet parameter
     //  YAe  ...Eccentricity dependet parameter
@@ -135,36 +149,34 @@ inline bool spheroidGeometry2(double radius, double aspectRatio,  //inputs
     //  XCe  ...Eccentricity dependet parameter
     //  YCe  ...Eccentricity dependet parameter
     //  YHe  ...Eccentricity dependet parameter
-  
 
-    double ei(0.0), Le(0.0);  
-    bool   result = 
-          spheroidGeometry(radius, aspectRatio,  //inputs
-                                 ai,     bi,              //outputs
-                                 ei,    Le              //outputs
-                              );
+    double ei(0.0), Le(0.0);
+    bool   result = spheroidGeometry(radius, aspectRatio, //inputs
+                                     ai,     bi,          //outputs
+                                     ei,     Le           //outputs
+                                    );
     if(!result)
         return false;
-        
+
     XAe= 2.6666666666666666666666667
          *ei*ei*ei
-         /(-2.0*ei+(1.0+ei*ei)*Le); 
+         /(-2.0*ei+(1.0+ei*ei)*Le);
     YAe= 5.333333333333333333333333333
         *ei*ei*ei
-        /(2.0*ei+(3*ei*ei-1.0)*Le);
+        /(2.0*ei+(3.0*ei*ei-1.0)*Le);
     XCe= 1.333333333333333333333333333
-        *ei*ei*ei
-        *(1.0-ei*ei)
-        /(2.0*ei-(1.0-ei*ei)*Le);
+         *ei*ei*ei
+         *(1.0-ei*ei)
+         /(2.0*ei-(1.0-ei*ei)*Le);
     YCe= 1.3333333333333333333333
-        *ei*ei*ei
-        *(2.0-ei*ei)
-        /(-2.0*ei+(1.0+ei*ei)*Le);
+         *ei*ei*ei
+         *(2.0-ei*ei)
+         /(-2.0*ei+(1.0+ei*ei)*Le);
     YHe= 1.3333333333333333333333
-        *ei*ei*ei*ei*ei
-        /(-2.0*ei+(1.0+ei*ei)*Le);
+         *ei*ei*ei*ei*ei
+         /(-2.0*ei+(1.0+ei*ei)*Le);
 
-	return true; 		    
+    return true;
 
 }
 
@@ -201,28 +213,28 @@ inline void multiply333(double scalar, double tensor[3][3][3] )
 inline void permutationTensor(double tensor[3][3][3] )
 {
     zeroize333(tensor);
-    tensor[0][1][2] = 1.0; 
-    tensor[1][2][0] = 1.0; 
-    tensor[2][0][1] = 1.0; 
-    tensor[0][2][1] =-1.0; 
-    tensor[2][1][0] =-1.0; 
-    tensor[1][0][2] =-1.0; 
+    tensor[0][1][2] = 1.0;
+    tensor[1][2][0] = 1.0;
+    tensor[2][0][1] = 1.0;
+    tensor[0][2][1] =-1.0;
+    tensor[2][1][0] =-1.0;
+    tensor[1][0][2] =-1.0;
 }
 
 
 
 //--------------------------------------------------------------------
-// Compute a dot product of the permutation tensor and 
+// Compute a dot product of the permutation tensor and
 // then a dyadic product of with a vector
-inline bool permutationDotDyadic( 
-                                 double vector[3], 
-                                 double  tensor[3][3][3]      
+inline bool permutationDotDyadic(
+                                 double vector[3],
+                                 double  tensor[3][3][3]
                                 )
 {
     //Generate permutation tensor
     double permutationT[3][3][3];
     permutationTensor(permutationT);
-    
+
     //Step 1: compute dot prodcut of permutation tensor
     double permutationDotProd[3][3];
     zeroize33(permutationDotProd);
@@ -237,20 +249,20 @@ inline bool permutationDotDyadic(
         for(int iY=0; iY<3; iY++)
             for(int iZ=0; iZ<3; iZ++)
                 tensor[iX][iY][iZ] = permutationDotProd[iX][iY]
-                                   * vector[iZ]; 
+                                   * vector[iZ];
     return true;
 
 }
 
 //--------------------------------------------------------------------
 // Compute a dot and dyadic product of with a vector
-inline bool doubleDotTensor333Tensor33(double tensor333[3][3][3], 
+inline bool doubleDotTensor333Tensor33(double tensor333[3][3][3],
                                        double tensor33[3][3],
                                        double result[3]
                                       )
 {
     result[0]=0.0;result[1]=0.0;result[2]=0.0;
-    
+
     for(int iX=0; iX<3; iX++)
         for(int iY=0; iY<3; iY++)
             for(int iZ=0; iZ<3; iZ++)
@@ -260,6 +272,6 @@ inline bool doubleDotTensor333Tensor33(double tensor333[3][3][3],
 }
 
 
-}; //end of namespace
+} //end of namespace
 
 #endif

src/lagrangian/cfdemParticle/cfdTools/setupProbeModel.H

@@ -1,3 +1,2 @@
- //set probeModel parameters for this force model
- particleCloud_.probeM().setOutputFile();
- particleCloud_.probeM().setCounter();
+//set probeModel parameters for this force model
+if (probeIt_) { particleCloud_.probeM().setOutputFile(typeName+".logDat"); }

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -5,7 +5,8 @@
     www.cfdem.com
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
-                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright 2012-2015 DCS Computing GmbH, Linz
+                                Copyright 2015-     JKU Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -33,12 +34,12 @@ Description
 #ifndef versionInfo_H
 #define versionInfo_H
 
-word CFDEMversion="cfdem-2.8.2";
-word compatibleLIGGGHTSversion="3.0.6";
-word OFversion="2.3.x-commit-4d6f4a3115ff76ec4154c580eb041bc95ba4ec09";
+word CFDEMversion="PFM 17.08";
+word compatibleLIGGGHTSversion="PFM 17.08";
+word OFversion="4.x";
 
 Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;
 Info << ", compatible to LIGGGHTS version: " << compatibleLIGGGHTSversion <<  endl;
-Info << ", compatible to OF version and build: " << OFversion <<  endl;
+Info << ", compatible to OpenFOAM version: " << OFversion <<  endl;
 
 #endif

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -44,9 +44,13 @@ Description
 #include "clockModel.H"
 #include "smoothingModel.H"
 #include "liggghtsCommandModel.H"
+#include "otherForceModel.H"
+
+namespace Foam
+{
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-Foam::cfdemCloud::cfdemCloud
+cfdemCloud::cfdemCloud
 (
     const fvMesh& mesh
 )
@@ -77,6 +81,8 @@ Foam::cfdemCloud::cfdemCloud
     solveFlow_(true),
     verbose_(false),
     ignore_(false),
+    allowCFDsubTimestep_(true),
+    limitDEMForces_(false),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
     velocities_(NULL),
@@ -93,11 +99,13 @@ Foam::cfdemCloud::cfdemCloud
     particleVolumes_(NULL),
     particleV_(NULL),
     numberOfParticles_(0),
+    d32_(-1),
     numberOfParticlesChanged_(false),
     arraysReallocated_(false),
     forceModels_(couplingProperties_.lookup("forceModels")),
     momCoupleModels_(couplingProperties_.lookup("momCoupleModels")),
     liggghtsCommandModelList_(liggghtsCommandDict_.lookup("liggghtsCommandModels")),
+    otherForceModels_(couplingProperties_.lookupOrDefault<wordList>("otherForceModels",wordList(0))),
     turbulenceModelType_(couplingProperties_.lookup("turbulenceModelType")),
     cg_(1.),
     cgOK_(true),
@@ -107,7 +115,7 @@ Foam::cfdemCloud::cfdemCloud
     treatVoidCellsAsExplicitForce_(false),
     useDDTvoidfraction_(false),
     ddtVoidfraction_
-    (   
+    (
         IOobject
         (
             "ddtVoidfraction",
@@ -119,17 +127,10 @@ Foam::cfdemCloud::cfdemCloud
         mesh,
         dimensionedScalar("zero", dimensionSet(0,0,-1,0,0), 0)  // 1/s
     ),
+    checkPeriodicCells_(false),
     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_
         )
@@ -197,7 +198,7 @@ Foam::cfdemCloud::cfdemCloud
         clockModel::New
         (
             couplingProperties_,
-            *this
+            mesh.time()
         )
     ),
     smoothingModel_
@@ -226,29 +227,31 @@ Foam::cfdemCloud::cfdemCloud
         solveFlow_=Switch(couplingProperties_.lookup("solveFlow"));
     if (couplingProperties_.found("imExSplitFactor"))
         imExSplitFactor_ = readScalar(couplingProperties_.lookup("imExSplitFactor"));
+
+    if(imExSplitFactor_ > 1.0)
+        FatalError << "You have set imExSplitFactor > 1 in your couplingProperties. Must be <= 1."
+                   << abort(FatalError);
+    if(imExSplitFactor_ < 0.0)
+        FatalError << "You have set imExSplitFactor < 0 in your couplingProperties. Must be >= 0."
+                   << abort(FatalError);
+
     if (couplingProperties_.found("treatVoidCellsAsExplicitForce"))
         treatVoidCellsAsExplicitForce_ = readBool(couplingProperties_.lookup("treatVoidCellsAsExplicitForce"));
     if (couplingProperties_.found("verbose")) verbose_=true;
     if (couplingProperties_.found("ignore")) ignore_=true;
+    if (couplingProperties_.found("limitDEMForces"))
+    {
+        limitDEMForces_=true;
+        maxDEMForce_ = readScalar(couplingProperties_.lookup("limitDEMForces"));
+    }
     if (turbulenceModelType_=="LESProperties")
         Info << "WARNING - LES functionality not yet tested!" << endl;
 
     if (couplingProperties_.found("useDDTvoidfraction"))
         useDDTvoidfraction_=true;
-    else        
+    else
         Info << "ignoring ddt(voidfraction)" << endl;
 
-    forceModel_ = new autoPtr<forceModel>[nrForceModels()];
-    for (int i=0;i<nrForceModels();i++)
-    {
-        forceModel_[i] = forceModel::New
-        (
-            couplingProperties_,
-            *this,
-            forceModels_[i]
-        );
-    }
-
     momCoupleModel_ = new autoPtr<momCoupleModel>[momCoupleModels_.size()];
     for (int i=0;i<momCoupleModels_.size();i++)
     {
@@ -260,6 +263,17 @@ Foam::cfdemCloud::cfdemCloud
         );
     }
 
+    forceModel_ = new autoPtr<forceModel>[nrForceModels()];
+    for (int i=0;i<nrForceModels();i++)
+    {
+        forceModel_[i] = forceModel::New
+        (
+            couplingProperties_,
+            *this,
+            forceModels_[i]
+        );
+    }
+
     // run liggghts commands from cfdem
     liggghtsCommand_ = new autoPtr<liggghtsCommandModel>[liggghtsCommandModelList_.size()];
     for (int i=0;i<liggghtsCommandModelList_.size();i++)
@@ -273,12 +287,60 @@ Foam::cfdemCloud::cfdemCloud
         );
     }
 
+    otherForceModel_ = new autoPtr<otherForceModel>[otherForceModels_.size()];
+    for (int i=0;i<otherForceModels_.size();i++)
+    {
+        otherForceModel_[i] = otherForceModel::New
+        (
+            couplingProperties_,
+            *this,
+            otherForceModels_[i]
+        );
+    }
+
     dataExchangeM().setCG();
-    if (!cgOK_ && cg_ > 1) FatalError<< "at least one of your models is not fit for cg !!!"<< abort(FatalError); 
+    Switch cgWarnOnly_(couplingProperties_.lookupOrDefault<Switch>("cgWarnOnly", true));
+    if (!cgOK_ && cg_ > 1)
+    {
+        if (cgWarnOnly_)
+            Warning << "at least one of your models is not fit for cg !!!" << endl;
+        else
+            FatalError << "at least one of your models is not fit for cg !!!" << abort(FatalError);
+    }
+
+    // check if simulation is a fully periodic box
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    int nPatchesCyclic(0);
+    int nPatchesNonCyclic(0);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        if (isA<cyclicPolyPatch>(pp) || isA<cyclicAMIPolyPatch>(pp))
+            ++nPatchesCyclic;
+        else if (!isA<processorPolyPatch>(pp))
+            ++nPatchesNonCyclic;
+    }
+
+    if (nPatchesNonCyclic == 0)
+    {
+        checkPeriodicCells_ = true;
+    }
+
+    //hard set checkperiodic cells if wished
+    if(this->couplingProperties().found("checkPeriodicCells"))
+    {
+        checkPeriodicCells_ = couplingProperties().lookupOrDefault<Switch>("checkPeriodicCells", checkPeriodicCells_);
+    }
+
+    if (nPatchesCyclic > 0 && nPatchesNonCyclic > 0)
+    {
+        if (verbose_) Info << "nPatchesNonCyclic=" << nPatchesNonCyclic << ", nPatchesCyclic=" << nPatchesCyclic << endl;
+        Warning << "Periodic handing is disabled because the domain is not fully periodic!\n" << endl;
+    }
 }
 
 // * * * * * * * * * * * * * * * * Destructors  * * * * * * * * * * * * * * //
-Foam::cfdemCloud::~cfdemCloud()
+cfdemCloud::~cfdemCloud()
 {
     clockM().evalPar();
     clockM().normHist();
@@ -297,8 +359,9 @@ Foam::cfdemCloud::~cfdemCloud()
     dataExchangeM().destroy(particleVolumes_,1);
     dataExchangeM().destroy(particleV_,1);
 }
+
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void Foam::cfdemCloud::getDEMdata()
+void cfdemCloud::getDEMdata()
 {
     dataExchangeM().getData("radius","scalar-atom",radii_);
     dataExchangeM().getData("x","vector-atom",positions_);
@@ -308,7 +371,7 @@ void Foam::cfdemCloud::getDEMdata()
         dataExchangeM().getData("dragAcc","vector-atom",fAcc_); // array is used twice - might be necessary to clean it first
 }
 
-void Foam::cfdemCloud::giveDEMdata()
+void cfdemCloud::giveDEMdata()
 {
     if(forceM(0).coupleForce())
     {
@@ -328,7 +391,7 @@ void Foam::cfdemCloud::giveDEMdata()
 
 // * * * * * * * * * * * * * * * protected Member Functions  * * * * * * * * * * * * * //
 
-void Foam::cfdemCloud::setNumberOfParticles(int nP)
+void cfdemCloud::setNumberOfParticles(int nP)
 {
     if(nP != numberOfParticles())
     {
@@ -337,12 +400,12 @@ void Foam::cfdemCloud::setNumberOfParticles(int nP)
     }
 }
 
-void Foam::cfdemCloud::findCells()
+void cfdemCloud::findCells()
 {
     locateM().findCell(NULL,positions_,cellIDs_,numberOfParticles());
 }
 
-void Foam::cfdemCloud::setForces()
+void cfdemCloud::setForces()
 {
     resetArray(fluidVel_,numberOfParticles(),3);
     resetArray(impForces_,numberOfParticles(),3);
@@ -350,9 +413,22 @@ void Foam::cfdemCloud::setForces()
     resetArray(DEMForces_,numberOfParticles(),3);
     resetArray(Cds_,numberOfParticles(),1);
     for (int i=0;i<cfdemCloud::nrForceModels();i++) cfdemCloud::forceM(i).setForce();
+
+    if (limitDEMForces_)
+    {
+        scalar maxF = 0.0;
+        for (int index = 0;index <  numberOfParticles(); ++index)
+        {
+            scalar F = mag(expForce(index));
+            if (F > maxF) maxF = F;
+            if (F > maxDEMForce_)
+              for(int i=0;i<3;i++) DEMForces_[index][i] *= maxDEMForce_/F;
+        }
+        Info << "largest particle-fluid interaction on particle: " << maxF << endl;
+    }
 }
 
-void Foam::cfdemCloud::setParticleForceField()
+void cfdemCloud::setParticleForceField()
 {
     averagingM().setVectorSum
     (
@@ -370,7 +446,7 @@ void Foam::cfdemCloud::setParticleForceField()
     );
 }
 
-void Foam::cfdemCloud::setVectorAverages()
+void cfdemCloud::setVectorAverages()
 {
     if(verbose_) Info << "- setVectorAverage(Us,velocities_,weights_)" << endl;
     averagingM().setVectorAverage
@@ -383,13 +459,16 @@ void Foam::cfdemCloud::setVectorAverages()
     );
     if(verbose_) Info << "setVectorAverage done." << endl;
 }
+
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
-void Foam::cfdemCloud::checkCG(bool ok)
+
+void cfdemCloud::checkCG(bool ok)
 {
     if(!cgOK_) return;
     if(!ok) cgOK_ = ok;
 }
-void Foam::cfdemCloud::setPos(double**& pos)
+
+void cfdemCloud::setPos(double**& pos)
 {
     for(int index = 0;index <  numberOfParticles(); ++index)
     {
@@ -398,79 +477,71 @@ void Foam::cfdemCloud::setPos(double**& pos)
         }
     }
 }
+
 // * * * * * * * * * * * * * * * ACCESS  * * * * * * * * * * * * * //
 
-label Foam::cfdemCloud::particleCell(int index)
+label cfdemCloud::particleCell(int index) const
 {
-    label cellI = cellIDs()[index][0];
-    return cellI;
+    return cellIDs()[index][0];
 }
 
-vector Foam::cfdemCloud::position(int index)
+vector cfdemCloud::position(int index) const
 {
-    vector pos;
-    for(int i=0;i<3;i++) pos[i] = positions()[index][i];
-    return pos;
+    return vector(positions()[index][0],positions()[index][1],positions()[index][2]);
 }
 
-vector Foam::cfdemCloud::velocity(int index)
+vector cfdemCloud::velocity(int index) const
 {
-    vector vel;
-    for(int i=0;i<3;i++) vel[i] = velocities()[index][i];
-    return vel;
+    return vector(velocities()[index][0],velocities()[index][1],velocities()[index][2]);
 }
 
-vector Foam::cfdemCloud::fluidVel(int index)
+vector cfdemCloud::expForce(int index) const
 {
-    vector vel;
-    for(int i=0;i<3;i++) vel[i] = fluidVels()[index][i];
-    return vel;
+    return vector(DEMForces()[index][0],DEMForces()[index][1],DEMForces()[index][2]);
 }
 
-const forceModel& Foam::cfdemCloud::forceM(int i)
+vector cfdemCloud::fluidVel(int index) const
+{
+    return vector(fluidVels()[index][0],fluidVels()[index][1],fluidVels()[index][2]);
+}
+
+const forceModel& cfdemCloud::forceM(int i)
 {
     return forceModel_[i];
 }
 
-int Foam::cfdemCloud::nrForceModels()
+label cfdemCloud::nrForceModels() const
 {
     return forceModels_.size();
 }
 
-scalar Foam::cfdemCloud::voidfraction(int index)
+label cfdemCloud::nrMomCoupleModels() const
+{
+    return momCoupleModels_.size();
+}
+
+scalar cfdemCloud::voidfraction(int index) const
 {
     return voidfractions()[index][0];
 }
 
-label Foam::cfdemCloud::liggghtsCommandModelIndex(word name)
+label cfdemCloud::liggghtsCommandModelIndex(word name) const
 {
-    int index=-1;
     forAll(liggghtsCommandModelList_,i)
     {
         if(liggghtsCommand()[i]().name() == name)
         {
-            index = i;
-            break;
+            return i;
         }
     }
-    return index;
-}
-
-std::vector<double*>* Foam::cfdemCloud::getVprobe()
-{
- return probeModel_->getVprobe();
-}
-
-std::vector<double>* Foam::cfdemCloud::getSprobe()
-{
- return probeModel_->getSprobe();
+    return -1;
 }
 
 // * * * * * * * * * * * * * * * WRITE  * * * * * * * * * * * * * //
 
 // * * *   write cfdemCloud internal data   * * * //
 
-bool Foam::cfdemCloud::evolve
+bool cfdemCloud::evolve
 (
     volScalarField& alpha,
     volVectorField& Us,
@@ -480,20 +551,20 @@ bool Foam::cfdemCloud::evolve
     numberOfParticlesChanged_ = false;
     arraysReallocated_=false;
     bool doCouple=false;
-    probeModel_->clearProbes();
 
     if(!ignore())
     {
-        if (dataExchangeM().couple())
+        if (dataExchangeM().doCoupleNow())
         {
             Info << "\n Coupling..." << endl;
+            dataExchangeM().couple(0);
             doCouple=true;
 
             // reset vol Fields
             clockM().start(16,"resetVolFields");
             if(verbose_)
             {
-                Info << "couplingStep:" << dataExchangeM().couplingStep() 
+                Info << "couplingStep:" << dataExchangeM().couplingStep()
                      << "\n- resetVolFields()" << endl;
             }
             averagingM().resetVectorAverage(averagingM().UsPrev(),averagingM().UsNext(),false);
@@ -531,24 +602,28 @@ bool Foam::cfdemCloud::evolve
             setVectorAverages();
 
 
-            //Smoothen "next" fields            
+            //Smoothen "next" fields
             smoothingM().dSmoothing();
             smoothingM().smoothen(voidFractionM().voidFractionNext());
 
             //only smoothen if we use implicit force coupling in cells void of particles
-            //because we need unsmoothened Us field to detect cells for explicit 
+            //because we need unsmoothened Us field to detect cells for explicit
             //force coupling
             if(!treatVoidCellsAsExplicitForce())
                 smoothingM().smoothenReferenceField(averagingM().UsNext());
-            
+
             clockM().stop("setVectorAverage");
         }
-        
+
         //============================================
-        //CHECK JUST TIME-INTERPOATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD 
+        //CHECK JUST TIME-INTERPOATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD
         //      IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED
         //      QUANTITIES AT THE GRID!
         Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;
+        if(dataExchangeM().timeStepFraction() > 1.0000001)
+        {
+            FatalError << "cfdemCloud::dataExchangeM().timeStepFraction()>1: Do not do this, since dangerous. This might be due to the fact that you used a adjustable CFD time step. Please use a fixed CFD time step." << abort(FatalError);
+        }
         clockM().start(24,"interpolateEulerFields");
 
         // update voidFractionField
@@ -591,6 +666,8 @@ bool Foam::cfdemCloud::evolve
             clockM().start(23,"giveDEMdata");
             giveDEMdata();
             clockM().stop("giveDEMdata");
+
+            dataExchangeM().couple(1);
         }//end dataExchangeM().couple()
 
 
@@ -607,7 +684,7 @@ bool Foam::cfdemCloud::evolve
     return doCouple;
 }
 
-bool Foam::cfdemCloud::reAllocArrays() const
+bool cfdemCloud::reAllocArrays()
 {
     if(numberOfParticlesChanged_ && !arraysReallocated_)
     {
@@ -622,7 +699,7 @@ bool Foam::cfdemCloud::reAllocArrays() const
         dataExchangeM().allocateArray(Cds_,0.,1);
         dataExchangeM().allocateArray(radii_,0.,1);
         dataExchangeM().allocateArray(voidfractions_,1.,voidFractionM().maxCellsPerParticle());
-        dataExchangeM().allocateArray(cellIDs_,-1.,voidFractionM().maxCellsPerParticle());
+        dataExchangeM().allocateArray(cellIDs_,-1,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleV_,0.,1);
@@ -632,29 +709,6 @@ bool Foam::cfdemCloud::reAllocArrays() const
     return false;
 }
 
-bool Foam::cfdemCloud::reAllocArrays(int nP, bool forceRealloc) const
-{
-    if( (numberOfParticlesChanged_ && !arraysReallocated_) || forceRealloc)
-    {
-        // get arrays of new length
-        dataExchangeM().allocateArray(positions_,0.,3,nP);
-        dataExchangeM().allocateArray(velocities_,0.,3,nP);
-        dataExchangeM().allocateArray(fluidVel_,0.,3,nP);
-        dataExchangeM().allocateArray(impForces_,0.,3,nP);
-        dataExchangeM().allocateArray(expForces_,0.,3,nP);
-        dataExchangeM().allocateArray(DEMForces_,0.,3,nP);
-        dataExchangeM().allocateArray(Cds_,0.,1,nP);
-        dataExchangeM().allocateArray(radii_,0.,1,nP);
-        dataExchangeM().allocateArray(voidfractions_,1.,voidFractionM().maxCellsPerParticle(),nP);
-        dataExchangeM().allocateArray(cellIDs_,0.,voidFractionM().maxCellsPerParticle(),nP);
-        dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle(),nP);
-        dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle(),nP);
-        arraysReallocated_ = true;
-        return true;
-    }
-    return false;
-}
-
 tmp<fvVectorMatrix> cfdemCloud::divVoidfractionTau(volVectorField& U,volScalarField& voidfraction) const
 {
     return
@@ -687,7 +741,7 @@ void cfdemCloud::calcDdtVoidfraction(volScalarField& voidfraction) const
 /*tmp<fvVectorMatrix> cfdemCloud::ddtVoidfractionU(volVectorField& U,volScalarField& voidfraction) const
 {
     if (dataExchangeM().couplingStep() <= 2) return fvm::ddt(U);
-    
+
     return fvm::ddt(voidfraction,U);
 }*/
 
@@ -725,8 +779,17 @@ void cfdemCloud::resetArray(double**& array,int length,int width,double resetVal
         }
     }
 }
-// * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
-#include "cfdemCloudIO.C"
+void cfdemCloud::otherForces(volVectorField& forcefield)
+{
+  forcefield.primitiveFieldRef() = vector::zero;
+  forcefield.boundaryFieldRef() = vector::zero;
+  for (int i=0;i<otherForceModels_.size();i++)
+      forcefield += otherForceModel_[i]().exportForceField();
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
 
 // ************************************************************************* //

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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,6 +69,7 @@ class smoothingModel;
 class momCoupleModel;
 class meshMotionModel;
 class liggghtsCommandModel;
+class otherForceModel;
 
 /*---------------------------------------------------------------------------*\
                            Class cfdemCloud Declaration
@@ -95,6 +92,12 @@ protected:
 
     bool ignore_;
 
+    bool allowCFDsubTimestep_;
+
+    bool limitDEMForces_;
+
+    scalar maxDEMForce_;
+
     const word modelType_;
 
     mutable double **positions_;
@@ -117,7 +120,7 @@ protected:
 
     mutable double **voidfractions_;
 
-    mutable double **cellIDs_;
+    mutable int **cellIDs_;
 
     mutable double **particleWeights_;
 
@@ -127,6 +130,8 @@ protected:
 
     int numberOfParticles_;
 
+    scalar d32_;
+
     bool numberOfParticlesChanged_;
 
     mutable bool arraysReallocated_;
@@ -137,6 +142,8 @@ protected:
 
     const wordList liggghtsCommandModelList_;
 
+    const wordList otherForceModels_;
+
     const word turbulenceModelType_;
 
     mutable scalar cg_;
@@ -155,15 +162,9 @@ 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
+    mutable Switch checkPeriodicCells_;
+
+    const turbulenceModel& turbulence_;
 
     autoPtr<forceModel>* forceModel_;
 
@@ -189,6 +190,8 @@ protected:
 
     autoPtr<liggghtsCommandModel>* liggghtsCommand_;
 
+    autoPtr<otherForceModel>* otherForceModel_;
+
 // Protected member functions
     virtual void getDEMdata();
 
@@ -225,41 +228,49 @@ public:
 // public Member Functions
 
     // Access
+       bool allowCFDsubTimestep() { return allowCFDsubTimestep_; }
+
+       void setAllowCFDsubTimestep(bool b) { allowCFDsubTimestep_ = b; }
+
        void checkCG(bool);
 
        void setPos(double **&);
 
-       word modelType(){ return modelType_; };
+       const word& modelType() const { return modelType_; }
 
-       label particleCell(int);
+       label particleCell(int) const;
 
-       vector position(int);
+       vector position(int) const;
 
-       vector velocity(int);
+       vector velocity(int) const;
 
-       vector fluidVel(int);
+       vector expForce(int) const;
+
+       vector fluidVel(int) const;
 
        virtual const forceModel& forceM(int);
 
-       virtual int nrForceModels();
+       virtual label nrForceModels() const;
 
-       scalar voidfraction(int);
+       virtual label nrMomCoupleModels() const;
 
-       label liggghtsCommandModelIndex(word);
+       scalar voidfraction(int) const;
 
-       inline void setCG(double) const;
+       label liggghtsCommandModelIndex(word) const;
 
-       inline const scalar& cg() const;
+       inline void setCG(double);
 
-       inline const bool& impDEMdrag() const;
+       inline scalar cg() const;
 
-       inline const bool& impDEMdragAcc() const;
+       inline bool impDEMdrag() const;
 
-       inline const scalar& imExSplitFactor() const;
+       inline bool impDEMdragAcc() const;
 
-       inline const bool& treatVoidCellsAsExplicitForce() const;
+       inline scalar imExSplitFactor() const;
 
-       inline const bool& ignore() const;
+       inline bool treatVoidCellsAsExplicitForce() const;
+
+       inline bool ignore() const;
 
        inline const fvMesh& mesh() const;
 
@@ -291,38 +302,39 @@ public:
 
        inline void get_radii(double**&) const;
 
-       inline double ** cellIDs() const;
+       inline int ** cellIDs() const;
 
-       inline void get_cellIDs(double**&) const;
+       inline void get_cellIDs(int**&) const;
 
        inline double ** particleWeights() const;
 
-       virtual inline label body(int);
+       virtual inline label body(int) const;
 
-       virtual inline double particleVolume(int);
+       virtual inline double particleVolume(int) const;
 
-       inline scalar radius(int);
+       inline scalar radius(int) const;
 
-       virtual inline double d(int);
+       virtual inline double d(int) const;
 
+       inline scalar d32(bool recalc=true);
        virtual inline double dMin() {return -1;}
        virtual inline double dMax() {return -1;}
        virtual inline int minType() {return -1;}
        virtual inline int maxType() {return -1;}
        virtual inline bool multipleTypesDMax() {return false;}
        virtual inline bool multipleTypesDMin() {return false;}
-       virtual inline double ** particleDensity() const {return NULL;};
-       virtual inline int ** particleTypes() const {return NULL;};
-       virtual label particleType(label index) const {return -1;};
+       virtual inline double ** particleDensity() const {return NULL;}
+       virtual inline int ** particleTypes() const {return NULL;}
+       virtual label particleType(label index) const {return -1;}
 
        //access to the particle's rotation and torque data
-       virtual inline double ** DEMTorques() const {return NULL;};
-       virtual inline double ** omegaArray() const {return NULL;};
-       virtual vector omega(int) const {return Foam::vector(0,0,0);};
+       virtual inline double ** DEMTorques() const {return NULL;}
+       virtual inline double ** omegaArray() const {return NULL;}
+       virtual vector omega(int) const {return vector(0,0,0);}
 
        //access to the particles' orientation information
-       virtual inline double ** exArray() const {return NULL;};
-       virtual vector ex(int) const {return Foam::vector(0,0,0);};
+       virtual inline double ** exArray() const {return NULL;}
+       virtual vector ex(int) const {return vector(0,0,0);}
 
        //Detector if SRF module is enable or not
        virtual inline bool SRFOn(){return false;}
@@ -335,7 +347,7 @@ public:
 
        inline bool arraysReallocated() const;
 
-       inline const wordList& forceModels();
+       inline const wordList& forceModels() const;
 
        inline const voidFractionModel& voidFractionM() const;
 
@@ -343,11 +355,11 @@ public:
 
        inline const momCoupleModel& momCoupleM(int) const;
 
-       inline const dataExchangeModel& dataExchangeM() const;
+       inline dataExchangeModel& dataExchangeM();
 
        inline const IOModel& IOM() const;
 
-       inline const probeModel& probeM() const;
+       inline probeModel& probeM();
 
        inline const averagingModel& averagingM() const;
 
@@ -361,30 +373,21 @@ 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
         virtual bool evolve(volScalarField&,volVectorField&,volVectorField&);
 
-        virtual bool reAllocArrays() const;
-
-        virtual bool reAllocArrays(int nP, bool forceRealloc) const;  //force number of particles during reallocation
+        virtual void postFlow() {}
 
+        virtual bool reAllocArrays();
 
     // IO
-        void writeScalarFieldToTerminal(double**&);
+        void writeScalarFieldToTerminal(double**&) const;
 
-        void writeVectorFieldToTerminal(double**&);
+        void writeVectorFieldToTerminal(double**&) const;
 
     // functions
         tmp<fvVectorMatrix> divVoidfractionTau(volVectorField& ,volScalarField&) const;
@@ -399,9 +402,9 @@ public:
 
         void resetArray(double**&,int,int,double resetVal=0.);
 
-        std::vector<double*>* getVprobe();
-         
-        std::vector<double>* getSprobe();
+        void otherForces(volVectorField&);
+
+        bool checkPeriodicCells() { return checkPeriodicCells_; }
 };
 
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -44,38 +44,38 @@ namespace Foam
 {
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-inline void cfdemCloud::setCG(double cg) const
-{ 
+inline void cfdemCloud::setCG(double cg)
+{
     cg_ = cg;
     Info << "cg is set to: " << cg_ << endl;
-};
+}
 
-inline const bool& cfdemCloud::impDEMdrag() const
-{ 
-    return impDEMdrag_;
-};
-
-inline const bool& cfdemCloud::impDEMdragAcc() const
-{ 
-    return impDEMdragAcc_;
-};
-
-inline const scalar& cfdemCloud::imExSplitFactor() const
+inline bool cfdemCloud::impDEMdrag() const
 {
-	    return imExSplitFactor_;
-};
+    return impDEMdrag_;
+}
 
-inline const bool& cfdemCloud::treatVoidCellsAsExplicitForce() const
+inline bool cfdemCloud::impDEMdragAcc() const
+{
+    return impDEMdragAcc_;
+}
+
+inline scalar cfdemCloud::imExSplitFactor() const
+{
+    return imExSplitFactor_;
+}
+
+inline bool cfdemCloud::treatVoidCellsAsExplicitForce() const
 {
     return treatVoidCellsAsExplicitForce_;
 }
 
-inline const scalar& cfdemCloud::cg() const
-{ 
+inline scalar cfdemCloud::cg() const
+{
     return cg_;
-};
+}
 
-inline const bool& cfdemCloud::ignore() const
+inline bool cfdemCloud::ignore() const
 {
     return ignore_;
 }
@@ -159,12 +159,12 @@ inline void cfdemCloud::get_radii(double **& values) const
     values=radii_;
 }
 
-inline double ** cfdemCloud::cellIDs() const
+inline int ** cfdemCloud::cellIDs() const
 {
     return cellIDs_;
 }
 
-inline void cfdemCloud::get_cellIDs(double **& values) const
+inline void cfdemCloud::get_cellIDs(int **& values) const
 {
 //    // make a copy of the array entries
 //    for (int i=0;i<numberOfParticles_;i++)
@@ -177,24 +177,43 @@ inline double ** cfdemCloud::particleWeights() const
     return particleWeights_;
 }
 
-inline label Foam::cfdemCloud::body(int index)
+inline label cfdemCloud::body(int index) const
 {
     return index;
 }
 
-inline double cfdemCloud::particleVolume(int index)
+inline double cfdemCloud::particleVolume(int index) const
 {
     return particleV_[index][0];
 }
 
-inline scalar cfdemCloud::radius(int index)
+inline scalar cfdemCloud::radius(int index) const
 {
     return radii_[index][0];
 }
 
-inline double cfdemCloud::d(int index)
+inline double cfdemCloud::d(int index) const
 {
-    return 2*radii_[index][0];
+    return 2.*radii_[index][0];
+}
+
+inline double cfdemCloud::d32(bool recalc)
+{
+    if(d32_<0 || recalc)
+    {
+        scalar Ntot(0);
+        scalar Dtot(0);
+        scalar r(0);
+        for(int index = 0;index <  numberOfParticles(); ++index)
+        {
+            r=radii_[index][0];
+            Ntot+=2*r*r*r;
+            Dtot+=r*r;
+        }
+        d32_=Ntot/Dtot;
+    }
+
+    return d32_;
 }
 
 inline int cfdemCloud::numberOfParticles() const
@@ -218,7 +237,7 @@ inline bool cfdemCloud::arraysReallocated() const
     return arraysReallocated_;
 }
 
-inline const wordList& cfdemCloud::forceModels()
+inline const wordList& cfdemCloud::forceModels() const
 {
     return forceModels_;
 }
@@ -233,9 +252,9 @@ inline const momCoupleModel& cfdemCloud::momCoupleM(int i) const
     return momCoupleModel_[i];
 }
 
-inline const dataExchangeModel& cfdemCloud::dataExchangeM() const
+inline dataExchangeModel& cfdemCloud::dataExchangeM()
 {
-    return dataExchangeModel_;
+    return dataExchangeModel_();
 }
 
 inline const IOModel& cfdemCloud::IOM() const
@@ -243,9 +262,9 @@ inline const IOModel& cfdemCloud::IOM() const
     return IOModel_;
 }
 
-inline const probeModel& cfdemCloud::probeM() const
+inline probeModel& cfdemCloud::probeM()
 {
-    return probeModel_;
+    return probeModel_();
 }
 
 inline const voidFractionModel& cfdemCloud::voidFractionM() const
@@ -283,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/cfdemCloud/cfdemCloudIO.C

@@ -32,6 +32,9 @@ Description
 
 #include "cfdemCloud.H"
 
+namespace Foam
+{
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 
@@ -45,7 +48,7 @@ Description
 
 // * * * * * * * * * * * * * * * IO * * * * * * * * * * * * * //
 
-void Foam::cfdemCloud::writeScalarFieldToTerminal(double**& array)
+void cfdemCloud::writeScalarFieldToTerminal(double**& array) const
 {
     // init double array
     for (int i=0; i<numberOfParticles(); i++)
@@ -54,7 +57,7 @@ void Foam::cfdemCloud::writeScalarFieldToTerminal(double**& array)
     }
 }
 
-void Foam::cfdemCloud::writeVectorFieldToTerminal(double**& array)
+void cfdemCloud::writeVectorFieldToTerminal(double**& array) const
 {
     // init double array
     for (int i=0; i<numberOfParticles(); i++)
@@ -65,4 +68,8 @@ void Foam::cfdemCloud::writeVectorFieldToTerminal(double**& array)
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 // ************************************************************************* //

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

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "cfdemCloudEnergy.H"
+#include "energyModel.H"
+#include "thermCondModel.H"
+#include "chemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+cfdemCloudEnergy::cfdemCloudEnergy
+(
+    const fvMesh& mesh
+)
+:
+    cfdemCloud(mesh),
+    energyModels_(couplingProperties_.lookup("energyModels")),
+    implicitEnergyModel_(false),
+    thermCondModel_
+    (
+        thermCondModel::New
+        (
+            couplingProperties_,
+            *this
+        )
+    ),
+    chemistryModel_
+    (
+        chemistryModel::New
+        (
+            couplingProperties_,
+            *this
+        )
+    )
+{
+    energyModel_ = new autoPtr<energyModel>[nrEnergyModels()];
+    for (int i=0;i<nrEnergyModels();i++)
+    {
+        energyModel_[i] = energyModel::New
+        (
+            couplingProperties_,
+            *this,
+            energyModels_[i]
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+cfdemCloudEnergy::~cfdemCloudEnergy()
+{
+
+}
+
+// * * * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
+
+void cfdemCloudEnergy::calcEnergyContributions()
+{
+    for (int i=0;i<nrEnergyModels();i++)
+        energyModel_[i]().calcEnergyContribution();
+}
+
+void cfdemCloudEnergy::speciesExecute()
+{
+        chemistryModel_().execute();
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+label cfdemCloudEnergy::nrEnergyModels() const
+{
+    return energyModels_.size();
+}
+
+bool cfdemCloudEnergy::implicitEnergyModel() const
+{
+    return implicitEnergyModel_;
+}
+
+const energyModel& cfdemCloudEnergy::energyM(int i)
+{
+    return energyModel_[i];
+}
+
+const chemistryModel& cfdemCloudEnergy::chemistryM()
+{
+    return chemistryModel_;
+}
+
+const thermCondModel& cfdemCloudEnergy::thermCondM()
+{
+    return thermCondModel_;
+}
+
+void cfdemCloudEnergy::energyContributions(volScalarField& Qsource)
+{
+    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.primitiveFieldRef()=0.0;
+    Qcoeff.boundaryFieldRef()=0.0;
+    for (int i=0;i<nrEnergyModels();i++)
+        energyM(i).addEnergyCoefficient(Qcoeff);
+}
+
+bool cfdemCloudEnergy::evolve
+(
+    volScalarField& alpha,
+    volVectorField& Us,
+    volVectorField& U
+)
+{
+    if (cfdemCloud::evolve(alpha, Us, U))
+    {
+        // calc energy contributions
+        // position 26 was already defined as Flow in clockModels and RhoPimpleChem solver.
+        clockM().start(27,"calcEnergyContributions");
+        if(verbose_) Info << "- calcEnergyContributions" << endl;
+        calcEnergyContributions();
+        if(verbose_) Info << "calcEnergyContributions done." << endl;
+        clockM().stop("calcEnergyContributions");
+
+        // execute chemical model species
+        clockM().start(28,"speciesExecute");
+        if(verbose_) Info << "- speciesExecute()" << endl;
+        speciesExecute();
+        if(verbose_) Info << "speciesExecute done" << endl;
+        clockM().stop("speciesExecute");
+
+	return true;
+    }
+    return false;
+}
+
+void cfdemCloudEnergy::postFlow()
+{
+    cfdemCloud::postFlow();
+    for (int i=0;i<nrEnergyModels();i++)
+        energyModel_[i]().postFlow();
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.H

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+
+Description
+    cloud class managing DEM data for CFD-DEM coupling with energy models,
+    e.g. for compressible, reacting flows
+
+Class
+    Foam::cfdemCloudEnergy derived from cfdemCloud
+
+SourceFiles
+    cfdemCloudEnergy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cfdemCloudEnergy_H
+#define cfdemCloudEnergy_H
+
+#include "cfdemCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+// forward declarations
+class energyModel;
+class thermCondModel;
+class chemistryModel;
+/*---------------------------------------------------------------------------*\
+                           Class cfdemCloudEnergy Declaration
+\*---------------------------------------------------------------------------*/
+
+class cfdemCloudEnergy
+:
+    public cfdemCloud
+{
+protected:
+
+    const wordList energyModels_;
+
+    bool implicitEnergyModel_;
+
+    const wordList chemistryModels_;
+
+    autoPtr<energyModel>* energyModel_;
+
+    autoPtr<thermCondModel> thermCondModel_;
+
+    autoPtr<chemistryModel> chemistryModel_;
+
+    void calcEnergyContributions();
+
+    void speciesExecute();
+
+public:
+
+    friend class energyModel;
+
+    // Constructors
+
+        //- Construct from components
+        cfdemCloudEnergy
+        (
+            const fvMesh& mesh
+        );
+
+    // Destructor
+
+        ~cfdemCloudEnergy();
+
+
+    // Member Functions
+
+
+       const energyModel& energyM(int);
+
+       const thermCondModel& thermCondM();
+
+       const chemistryModel& chemistryM();
+
+       label nrEnergyModels() const;
+
+       inline const wordList& energyModels() const;
+
+       bool implicitEnergyModel() const;
+
+       void energyContributions(volScalarField&);
+
+       void energyCoefficients(volScalarField&);
+
+       bool evolve(volScalarField&,volVectorField&,volVectorField&);
+
+       void postFlow();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "cfdemCloudEnergyI.H"
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergyI.H

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+inline const wordList& cfdemCloudEnergy::energyModels() const
+{
+    return energyModels_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C

@@ -36,7 +36,7 @@ Description
 #include "locateModel.H"
 #include "dataExchangeModel.H"
 #include "IOModel.H"
-#include "mpi.h"
+#include <mpi.h>
 #include "IOmanip.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -74,46 +74,49 @@ cfdemCloudIB::cfdemCloudIB
 
 cfdemCloudIB::~cfdemCloudIB()
 {
-    delete angularVelocities_;
+    dataExchangeM().destroy(angularVelocities_,3);
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-void Foam::cfdemCloudIB::getDEMdata()
+void cfdemCloudIB::getDEMdata()
 {
     cfdemCloud::getDEMdata();
     Info << "=== cfdemCloudIB::getDEMdata() === particle rotation not considered in CFD" << endl;
     //dataExchangeM().getData("omega","vector-atom",angularVelocities_);
 }
 
-bool Foam::cfdemCloudIB::reAllocArrays() const
+bool cfdemCloudIB::reAllocArrays()
 {
     if(cfdemCloud::reAllocArrays())
     {
-        dataExchangeM().allocateArray(angularVelocities_,0,3);
+        // get arrays of new length
+        dataExchangeM().allocateArray(angularVelocities_,0.,3);
+        return true;
     }
-    return true;
+    return false;
 }
 
-bool Foam::cfdemCloudIB::evolve()
+bool cfdemCloudIB::evolve()
 {
     numberOfParticlesChanged_ = false;
     arraysReallocated_=false;
     bool doCouple=false;
 
-    if (dataExchangeM().couple())
+    if (dataExchangeM().doCoupleNow())
     {
         Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;
+        dataExchangeM().couple(0);
         doCouple=true;
 
-//        Info << "skipLagrangeToEulerMapping_: " << skipLagrangeToEulerMapping_ 
+//        Info << "skipLagrangeToEulerMapping_: " << skipLagrangeToEulerMapping_
 //             << " haveEvolvedOnce_: " << haveEvolvedOnce_ << endl;
         if(!skipLagrangeToEulerMapping_ || !haveEvolvedOnce_)
         {
           if(verbose_) Info << "- getDEMdata()" << endl;
           getDEMdata();
           Info << "nr particles = " << numberOfParticles() << endl;
-        
+
           // search cellID of particles
           if(verbose_) Info << "- findCell()" << endl;
           locateM().findCell(NULL,positions_,cellIDs_,numberOfParticles());
@@ -124,7 +127,7 @@ bool Foam::cfdemCloudIB::evolve()
           voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_,particleV_);
           if(verbose_) Info << "setvoidFraction done." << endl;
         }
-        
+
         // set particles forces
         if(verbose_) Info << "- setForce(forces_)" << endl;
         for(int index = 0;index <  numberOfParticles_; ++index){
@@ -140,7 +143,9 @@ bool Foam::cfdemCloudIB::evolve()
         // write DEM data
         if(verbose_) Info << " -giveDEMdata()" << endl;
         giveDEMdata();
-        
+
+        dataExchangeM().couple(1);
+
         haveEvolvedOnce_=true;
     }
     Info << "evolve done." << endl;
@@ -153,7 +158,7 @@ bool Foam::cfdemCloudIB::evolve()
     return doCouple;
 }
 
-void Foam::cfdemCloudIB::calcVelocityCorrection
+void cfdemCloudIB::calcVelocityCorrection
 (
     volScalarField& p,
     volVectorField& U,
@@ -197,11 +202,11 @@ void Foam::cfdemCloudIB::calcVelocityCorrection
     (
         fvm::laplacian(phiIB) == fvc::div(U) + fvc::ddt(voidfraction)
     );
-     if(phiIB.needReference()) 
-     {
-         phiIBEqn.setReference(pRefCell_, pRefValue_);
-     }
-    
+    if(phiIB.needReference())
+    {
+        phiIBEqn.setReference(pRefCell_, pRefValue_);
+    }
+
     phiIBEqn.solve();
 
     U=U-fvc::grad(phiIB);
@@ -211,19 +216,17 @@ void Foam::cfdemCloudIB::calcVelocityCorrection
     p=p+phiIB/U.mesh().time().deltaT();  // do we have to  account for rho here?
     p.correctBoundaryConditions();
 
-     if (couplingProperties_.found("checkinterface"))
-       {
-          Info << "checking no-slip on interface..." << endl;
+    if (couplingProperties_.found("checkinterface"))
+    {
+        Info << "checking no-slip on interface..." << endl;
 //          #include "checkInterfaceVelocity.H" //TODO: check carefully!
-       }
+    }
 
 }
 
-vector Foam::cfdemCloudIB::angularVelocity(int index)
+vector cfdemCloudIB::angularVelocity(int index) const
 {
-    vector vel;
-    for(int i=0;i<3;i++) vel[i] = angularVelocities_[index][i]; 
-    return vel;
+    return vector(angularVelocities_[index][0],angularVelocities_[index][1],angularVelocities_[index][2]);
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H

@@ -61,7 +61,7 @@ protected:
     mutable double **angularVelocities_;
     label pRefCell_;
     scalar pRefValue_;
-    
+
     mutable bool haveEvolvedOnce_;
     mutable bool skipLagrangeToEulerMapping_;
 
@@ -83,19 +83,19 @@ public:
     // Member Functions
         void getDEMdata();
 
-        bool reAllocArrays() const;
-    
+        bool reAllocArrays();
+
         bool evolve();
 
         void calcVelocityCorrection(volScalarField&,volVectorField&,volScalarField&,volScalarField&); // this could be moved to an IB mom couple model
 
       // Access
-        vector angularVelocity(int);
+        vector angularVelocity(int) const;
 
         inline double ** angularVelocities() const
         {
             return angularVelocities_;
-        };
+        }
 
 };
 

src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.C

@@ -35,7 +35,7 @@ Description
 #include "locateModel.H"
 #include "dataExchangeModel.H"
 
-//#include "mpi.h" // only for debug reason
+//#include <mpi.h> // only for debug reason
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -105,13 +105,6 @@ cfdemCloudMS::~cfdemCloudMS()
     dataExchangeM().destroy(clumpVol_,1);
     dataExchangeM().destroy(clumpDH_,1);
     dataExchangeM().destroy(clumpWeights_,1);
-    //delete exCM_;
-    //delete eyCM_;
-    //delete ezCM_;
-    //delete SclumpCM_;
-    //delete scalingCM_;
-    //delete Cclump_ex_;
-    //delete Cclump_ey_;
     dataExchangeM().destroy(impForcesCM_,3);
     dataExchangeM().destroy(expForcesCM_,3);
     dataExchangeM().destroy(DEMForcesCM_,3);
@@ -135,14 +128,14 @@ void cfdemCloudMS::getDEMdata()
 
     //- save clump volume and mass
     double **typeDH(NULL);
-    dataExchangeM().allocateArray(typeDH,-1,1,nClumpTypes()+1);
+    dataExchangeM().allocateArray(typeDH,-1.,1,nClumpTypes()+1);
     if(manDHdev_) // use manually defined dH
     {
         for(int k = 1;k <= nClumpTypes(); k++)
             typeDH[k][0]=dHbyV_[k-1]*typeVol_[k];
     }
     else // calc dH from volAeqivalent shpere
-    {  
+    {
         for(int k = 1;k <= nClumpTypes(); k++)
             typeDH[k][0]=pow(typeVol_[k]*1.9099,1./3.); // 6/pi=1.9099 // calc a hydraulic diameter as d of vol equal sphere
     }
@@ -151,7 +144,7 @@ void cfdemCloudMS::getDEMdata()
     for(int ind = 0;ind < numberOfClumps(); ind++)
     {
         ct=clumpType()[0][ind];
-        clumpVol_[ind][0] = typeVol_[ct];      
+        clumpVol_[ind][0] = typeVol_[ct];
         clumpDH_[ind][0]=typeDH[ct][0];
         //Info << "ct=" << ct << endl;
         //Info << "clumpVol()[ind][0]=" << clumpVol()[ind][0] << endl;
@@ -169,7 +162,7 @@ void cfdemCloudMS::getDEMdata()
     //dataExchangeM().getScalarData("Cclump_ey",Cclump_ey_);   // cross section of the clump in ey normal direction
 }
 
-void Foam::cfdemCloudMS::giveDEMdata()
+void cfdemCloudMS::giveDEMdata()
 {
     /*for(int index = 0;index < numberOfClumps(); ++index)
     {
@@ -194,36 +187,29 @@ bool cfdemCloudMS::evolve
     return false;
 }
 
-bool cfdemCloudMS::reAllocArrays() const
+bool cfdemCloudMS::reAllocArrays()
 {
     if(cfdemCloud::reAllocArrays())
     {
         // get arrays of new length
-        dataExchangeM().allocateArray(positionsCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(velocitiesCM_,0,3,"nbodies");
+        dataExchangeM().allocateArray(positionsCM_,0.,3,"nbodies");
+        dataExchangeM().allocateArray(velocitiesCM_,0.,3,"nbodies");
         dataExchangeM().allocateArray(cellIDsCM_,-1,1,"nbodies");
         dataExchangeM().allocateArray(bodies_,0,1);
-        dataExchangeM().allocateArray(nrigids_,0,1,"nbodies");      
+        dataExchangeM().allocateArray(nrigids_,0,1,"nbodies");
         dataExchangeM().allocateArray(clumpType_,0,1,"nbodies");
-        dataExchangeM().allocateArray(clumpVol_,0,1,"nbodies");
+        dataExchangeM().allocateArray(clumpVol_,0.,1,"nbodies");
         dataExchangeM().allocateArray(clumpDH_,1.,1,"nbodies");
-        dataExchangeM().allocateArray(clumpWeights_,1,1,"nbodies");
-        //dataExchangeM().allocateArray(exCM_,0,3,"nbodies");
-        //dataExchangeM().allocateArray(eyCM_,0,3,"nbodies");
-        //dataExchangeM().allocateArray(ezCM_,0,3,"nbodies");
-        //dataExchangeM().allocateArray(SclumpCM_,0,3,nClumpTypes);
-        //dataExchangeM().allocateArray(scalingCM_,0,3,"nbodies");
-        //dataExchangeM().allocateArray(Cclump_ex_,0,3,nClumpTypes);
-        //dataExchangeM().allocateArray(Cclump_ey_,0,3,nClumpTypes);
-        dataExchangeM().allocateArray(impForcesCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(expForcesCM_,0,3,"nbodies");
-        dataExchangeM().allocateArray(DEMForcesCM_,0,3,"nbodies");
+        dataExchangeM().allocateArray(clumpWeights_,1.,1,"nbodies");
+        dataExchangeM().allocateArray(impForcesCM_,0.,3,"nbodies");
+        dataExchangeM().allocateArray(expForcesCM_,0.,3,"nbodies");
+        dataExchangeM().allocateArray(DEMForcesCM_,0.,3,"nbodies");
         return true;
     }
     return false;
 }
 
-void Foam::cfdemCloudMS::setNumberOfParticles(int nP)
+void cfdemCloudMS::setNumberOfParticles(int nP)
 {
     cfdemCloud::setNumberOfParticles(nP);
     int nC = dataExchangeM().getNumberOfClumps();
@@ -238,13 +224,13 @@ void Foam::cfdemCloudMS::setNumberOfParticles(int nP)
     numberOfClumps_ = min(numberOfParticles(),numberOfClumps_);
 }
 
-void Foam::cfdemCloudMS::findCells()
+void cfdemCloudMS::findCells()
 {
     cfdemCloud::findCells();
     locateM().findCell(NULL,positionsCM_,cellIDsCM_,numberOfClumps());
 }
 
-void Foam::cfdemCloudMS::setForces()
+void cfdemCloudMS::setForces()
 {
     resetArray(impForces_,numberOfParticles(),3);
     resetArray(expForces_,numberOfParticles(),3);
@@ -258,7 +244,7 @@ void Foam::cfdemCloudMS::setForces()
     for (int i=0;i<cfdemCloudMS::nrForceModels();i++) cfdemCloudMS::forceM(i).setForce();
 }
 
-void Foam::cfdemCloudMS::setParticleForceField()
+void cfdemCloudMS::setParticleForceField()
 {
     // set forces per particle
     cfdemCloud::setParticleForceField();
@@ -287,7 +273,7 @@ void Foam::cfdemCloudMS::setParticleForceField()
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 //             PUBLIC MEMBER FUNCTIONS
 
-const forceModel& Foam::cfdemCloudMS::forceM(int i)
+const forceModel& cfdemCloudMS::forceM(int i)
 {
     return forceModel_[i];
 }

src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.H

@@ -61,7 +61,7 @@ class cfdemCloudMS
 private:
     mutable double **positionsCM_;
     mutable double **velocitiesCM_;
-    mutable double **cellIDsCM_;
+    mutable int **cellIDsCM_;
     mutable int **bodies_;
     mutable int **nrigids_;
 
@@ -95,7 +95,6 @@ private:
 // Private member functions
     void getDEMdata();
     void giveDEMdata();
-    bool reAllocArrays() const;
     void setNumberOfParticles(int);
     void findCells();
     void setForces();
@@ -118,9 +117,9 @@ public:
 
     // Member Functions
       // Access
-       inline label body(int);
+       inline label body(int) const;
 
-       inline double particleVolume(int);
+       inline double particleVolume(int) const;
 
        inline vector positionCM(int);
 
@@ -128,9 +127,9 @@ public:
 
        inline label cellIDCM(int);
 
-       inline label nrigid(int);
+       inline label nrigid(int) const;
 
-       inline int nrForceModels();
+       inline label nrForceModels() const;
 
        inline double **& positionsCM() const;
 
@@ -164,6 +163,8 @@ public:
 
        const forceModel& forceM(int);
 
+       bool reAllocArrays();
+
        bool evolve(volScalarField&,volVectorField&,volVectorField&);
 };
 

src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMSI.H

@@ -29,17 +29,19 @@ Description
     and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 \*---------------------------------------------------------------------------*/
 
+#include "cfdemCloudMS.H"
+
 namespace Foam
 {
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label Foam::cfdemCloudMS::body(int index)
+inline label cfdemCloudMS::body(int index) const
 {
     return bodies_[0][index]-1;
 }
 
-inline double cfdemCloudMS::particleVolume(int index)
+inline double cfdemCloudMS::particleVolume(int index) const
 {
     int ind = body(index); // particle to clump ID
 
@@ -48,37 +50,37 @@ inline double cfdemCloudMS::particleVolume(int index)
         Vp=clumpVol_[ind][0];
     else
         Warning << "ind=" << ind << endl;
-    
+
     int nR(nrigid(ind));
     if(nR>0) Vp/=nR;
     return Vp;
 }
 
-inline vector Foam::cfdemCloudMS::positionCM(int index)
+inline vector cfdemCloudMS::positionCM(int index)
 {
     vector pos;
     for(int i=0;i<3;i++) pos[i] = positionsCM()[index][i];
     return pos;
 }
 
-inline vector Foam::cfdemCloudMS::velocityCM(int index)
+inline vector cfdemCloudMS::velocityCM(int index)
 {
     vector vel;
     for(int i=0;i<3;i++) vel[i] = velocitiesCM()[index][i];
     return vel;
 }
 
-inline label Foam::cfdemCloudMS::cellIDCM(int index)
+inline label cfdemCloudMS::cellIDCM(int index)
 {
     return cellIDsCM_[index][0];
 }
 
-inline label Foam::cfdemCloudMS::nrigid(int index)
+inline label cfdemCloudMS::nrigid(int index) const
 {
     return nrigids_[0][index];
 }
 
-inline int Foam::cfdemCloudMS::nrForceModels()
+inline label cfdemCloudMS::nrForceModels() const
 {
     return forceModels_.size();
 }

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

@@ -1,23 +0,0 @@
-#define version23  // currently being tested
-//#define version22  // currently being used
-//#define version21
-//#define version16ext
-//#define version15
-
-// choose comp/incomp
-//#define compre   // if comp is on - you must use Make/options_comp!
-
-//define multi sphere
-//#define multisphere
-
-// features of 2.1 work also in 2.3
-#if defined(version23)
-    #define version21
-    #define version221
-#endif
-
-// features of 2.1 work also in 2.2
-#if defined(version22)
-    #define version21
-    #define version221
-#endif

src/lagrangian/cfdemParticle/etc/bashrc

@@ -1,208 +0,0 @@
-#----------------------------------*-sh-*--------------------------------------
-# CFDEMcoupling
-# Christoph Goniva
-# June 2012
-#------------------------------------------------------------------------------
-#
-# Script
-#     etc/bashrc
-#
-# Description
-#     Startup file for cfdem exporting environment variables
-#     Sourced  ~/.bashrc
-#
-#- adapt this and add to ./bashrc 
-#- you can test the correctness using cfdemSystemTest.sh
-
-##================================================#
-##- source cfdem env vars
-#export CFDEM_VERSION=PUBLIC
-#export CFDEM_PROJECT_DIR=$HOME/CFDEM/CFDEMcoupling-$CFDEM_VERSION-$WM_PROJECT_VERSION
-#export CFDEM_SRC_DIR=$CFDEM_PROJECT_DIR/src
-#export CFDEM_SOLVER_DIR=$CFDEM_PROJECT_DIR/applications/solvers
-#export CFDEM_DOC_DIR=$CFDEM_PROJECT_DIR/doc
-#export CFDEM_UT_DIR=$CFDEM_PROJECT_DIR/applications/utilities
-#export CFDEM_TUT_DIR=$CFDEM_PROJECT_DIR/tutorials
-#export CFDEM_PROJECT_USER_DIR=$HOME/CFDEM/$LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION
-#export CFDEM_bashrc=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/bashrc
-#export CFDEM_LIGGGHTS_SRC_DIR=$HOME/LIGGGHTS/LIGGGHTS-PUBLIC/src
-#export CFDEM_LIGGGHTS_MAKEFILE_NAME=fedora_fpic
-#export CFDEM_LPP_DIR=$HOME/LIGGGHTS/mylpp/src
-#export CFDEM_PIZZA_DIR=$HOME/LIGGGHTS/PIZZA/gran_pizza_17Aug10/src
-#. $CFDEM_bashrc
-#================================================#
-#------------------------------------------------------------------------------
-
-#------------------------------------------------------------------------------
-#- export environment variables (adapt to your paths)
-#------------------------------------------------------------------------------
-
-#check if default lammps lib path should be used
-if [[ $CFDEM_LAMMPS_LIB_DIR == "" ]]; then
-    export CFDEM_LAMMPS_LIB_DIR=$CFDEM_LIGGGHTS_SRC_DIR/../lib/
-else
-    echo "using already defined CFDEM_LAMMPS_LIB_DIR=$CFDEM_LAMMPS_LIB_DIR ."
-fi
-
-#- LIGGGHTS lib name
-export CFDEM_LIGGGHTS_LIB_NAME=lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME
-
-#- CFDEM lib name
-export CFDEM_LIB_NAME=lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
-
-#- CFDEM compressible lib name
-export CFDEM_LIB_COMP_NAME=lagrangianCFDEMcomp-$CFDEM_VERSION-$WM_PROJECT_VERSION
-
-#-----------------------------------------------------
-# additional libraries
-
-#- LMP Many2Many lib path and makefile
-export CFDEM_Many2ManyLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
-export CFDEM_Many2ManyLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME
-
-#- LMP M2M lib path and makefile
-export CFDEM_M2MLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/library
-export CFDEM_M2MLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME
-
-#- LMP POEMS lib path and makefile
-export CFDEM_POEMSLIB_PATH=$CFDEM_LAMMPS_LIB_DIR/poems
-export CFDEM_POEMSLIB_MAKEFILENAME=g++
-
-#- LMP ASPHERE lib path and makefile
-export CFDEM_ASPHERELIB_PATH=$CFDEM_LAMMPS_LIB_DIR/poems
-export CFDEM_ASPHERELIB_MAKEFILENAME=g++
-
-#-C3PO library
-export C3PO_SRC_DIR=$CFDEM_SRC_DIR/c3po
-#-----------------------------------------------------
-
-#- path to test harness
-export CFDEM_TEST_HARNESS_PATH=$CFDEM_PROJECT_USER_DIR/log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
-
-#- path to libraries
-export CFDEM_LIB_DIR=$FOAM_USER_LIBBIN
-
-#- path to apps
-export CFDEM_APP_DIR=$FOAM_USER_APPBIN
-
-#- path to OF version flag file
-export CFDEM_OFVERSION_DIR=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/OFversion
-#------------------------------------------------------------------------------
-#- settings for lpp postproc tool
-#------------------------------------------------------------------------------
-
-#- nr of procs for lpp tool
-export CFDEM_LPP_NPROCS=4
-
-#- nr of procs for lpp tool
-export CFDEM_LPP_CHUNKSIZE=1
-
-#- shortcut to run lpp
-alias lpp='python -i $CFDEM_LPP_DIR/lpp.py --cpunum $CFDEM_LPP_NPROCS --chunksize $CFDEM_LPP_CHUNKSIZE'
-
-#------------------------------------------------------------------------------
-#- aliases for easy navigation (no changes necessary)
-#------------------------------------------------------------------------------
-#- shortcut to cfdem path
-alias cfdem='cd $CFDEM_PROJECT_DIR'
-
-#- shortcut to src path
-alias cfdemSrc='cd $CFDEM_SRC_DIR'
-
-#- shortcut to tutorial path
-alias cfdemTut='cd $CFDEM_TUT_DIR'
-
-#- shortcut to solver path
-alias cfdemSol='cd $CFDEM_SOLVER_DIR'
-
-#- shortcut to utilities path
-alias cfdemUt='cd $CFDEM_UT_DIR'
-
-#- shortcut to run path
-alias cfdemRun='cd $CFDEM_PROJECT_USER_DIR/run'
-
-#- shortcut to user solver path
-alias cfdemUsrSol='cd $CFDEM_PROJECT_USER_DIR/applications/solvers'
-
-#- shortcut to documentation path
-alias cfdemDoc='cd $CFDEM_DOC_DIR'
-
-#- shortcut to open the doxygen with firefox
-alias cfdemDox='firefox $CFDEM_DOC_DIR/doxygen/html/index.html'
-
-#- shortcut to LIGGGHTS path
-alias cfdemLIG='cd $CFDEM_LIGGGHTS_SRC_DIR'
-
-#- shortcut to system test
-alias cfdemSysTest='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cfdemSystemTest.sh'
-
-#- shortcut to pull LIGGGHTS
-alias cfdemPullLIG='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/pullLIGGGHTS.sh'
-
-#- shortcut to pull CFDEMcoupling
-alias cfdemPullCFDEMcoupling='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/pullCFDEMcoupling.sh'
-
-#- shortcut to clean CFDEM
-alias cfdemCleanCFDEM='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cleanCFDEMcoupling.sh'
-
-#- shortcut to compile LIGGGHTS + sublibraries
-alias cfdemCompLIG='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS.sh'
-
-#- shortcut to compile CFDEMcoupling +LIGGGHTS
-alias cfdemCompCFDEMall='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_all.sh'
-
-#- shortcut to compile CFDEMcoupling (src+solvers)
-alias cfdemCompCFDEM='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling.sh'
-
-#- shortcut to compile CFDEMcoupling src
-alias cfdemCompCFDEMsrc='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_src.sh'
-
-#- shortcut to compile CFDEMcoupling solvers
-alias cfdemCompCFDEMsol='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_sol.sh'
-
-#- shortcut to compile CFDEMcoupling utilities
-alias cfdemCompCFDEMuti='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh'
-
-#- shortcut to test basic tutorials
-alias cfdemTestTUT='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/testTutorials.sh'
-
-#- shortcut to visualize the clock model data
-alias vizClock='python $CFDEM_UT_DIR/vizClock/matPlot.py'
-
-#- recursive touch of current directory
-alias touchRec='find ./* -exec touch {} \;'
-
-#- shortcut to run liggghts in serial
-cfdemLiggghts() { $CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME < $1; }
-export -f cfdemLiggghts
-
-#- shortcut to run liggghts in parallel
-cfdemLiggghtsPar() { mpirun -np $2 $CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME < $1; }
-export -f cfdemLiggghtsPar
-
-#- shortcut to open files including a pattern
-cfdemGrep() { grep -rl "$1" ./* | xargs gedit; }
-export -f cfdemGrep
-
-#- shortcut lo list files in a directory
-#cfdemListFiles() { find $1 | sed s:""$1"":: > listOfFiles.txt; }   #leave out the dir iteslf in list
-cfdemListFiles() { find $1 > listOfFiles.txt; }                     #keep the dir in list
-export -f cfdemListFiles
-
-# check if the run directory exists
-if [ -d "$CFDEM_PROJECT_USER_DIR" ]; then
-    :
-else
-        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
-        read YN
-        if [ $YN == "y" ]; then
-            mkdir -p $CFDEM_PROJECT_USER_DIR
-            cd $CFDEM_PROJECT_USER_DIR
-            mkdir run
-            mkdir -p log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
-            mkdir -p applications/solvers
-        else
-            echo "aborted by user."
-            #exit
-        fi
-fi

src/lagrangian/cfdemParticle/etc/cshrc

@@ -1,185 +0,0 @@
-#----------------------------------*-sh-*--------------------------------------
-# CFDEMcoupling
-# Christoph Goniva
-# June 2012
-#------------------------------------------------------------------------------
-#
-# Script
-#     etc/cshrc
-#
-# Description
-#     Startup file for cfdem exporting environment variables
-#     Sourced  ~/.cshrc
-#
-#- adapt this and add to ./cshrc 
-#- you can test the correctness using cfdemSystemTest.sh
-
-##================================================#
-##- source cfdem env vars
-#setenv CFDEM_VERSION PUBLIC
-#setenv CFDEM_PROJECT_DIR $HOME/CFDEM/CFDEMcoupling-$CFDEM_VERSION-$WM_PROJECT_VERSION
-#setenv CFDEM_SRC_DIR $CFDEM_PROJECT_DIR/src/lagrangian/cfdemParticle
-#setenv CFDEM_SOLVER_DIR $CFDEM_PROJECT_DIR/applications/solvers
-#setenv CFDEM_DOC_DIR $CFDEM_PROJECT_DIR/doc
-#setenv CFDEM_UT_DIR $CFDEM_PROJECT_DIR/applications/utilities
-#setenv CFDEM_TUT_DIR $CFDEM_PROJECT_DIR/tutorials
-#setenv CFDEM_PROJECT_USER_DIR $HOME/CFDEM/$LOGNAME-$CFDEM_VERSION-$WM_PROJECT_VERSION
-#setenv CFDEM_bashrc $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cshrc
-#setenv CFDEM_LIGGGHTS_SRC_DIR $HOME/LIGGGHTS/LIGGGHTS-PUBLIC/src
-#setenv CFDEM_LIGGGHTS_MAKEFILE_NAME fedora_fpic
-#setenv CFDEM_LPP_DIR $HOME/LIGGGHTS/mylpp/src
-#setenv CFDEM_PIZZA_DIR $HOME/LIGGGHTS/PIZZA/gran_pizza_17Aug10/src
-#source $CFDEM_bashrc
-#================================================#
-#------------------------------------------------------------------------------
-
-#------------------------------------------------------------------------------
-#- export environment variables (adapt to your paths)
-#------------------------------------------------------------------------------
-
-#- LIGGGHTS lib name
-setenv CFDEM_LIGGGHTS_LIB_NAME lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME
-
-#- CFDEM lib name
-setenv CFDEM_LIB_NAME lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
-
-#- CFDEM compressible lib name
-setenv CFDEM_LIB_COMP_NAME lagrangianCFDEMcomp-$CFDEM_VERSION-$WM_PROJECT_VERSION
-
-#- LMP Many2Many lib path and makefile
-setenv CFDEM_Many2ManyLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
-setenv CFDEM_Many2ManyLIB_MAKEFILENAME $CFDEM_LIGGGHTS_MAKEFILE_NAME
-
-#- LMP M2M lib path
-setenv CFDEM_M2MLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/library
-setenv CFDEM_M2MLIB_MAKEFILENAME $CFDEM_LIGGGHTS_MAKEFILE_NAME
-
-#- LMP POEMS lib path
-setenv CFDEM_POEMSLIB_PATH $CFDEM_LIGGGHTS_SRC_DIR/../lib/poems
-setenv CFDEM_POEMSLIB_MAKEFILENAME g++
-
-#- path to test harness
-setenv CFDEM_TEST_HARNESS_PATH $CFDEM_PROJECT_USER_DIR/log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
-
-#- path to libraries
-setenv CFDEM_LIB_DIR $FOAM_USER_LIBBIN
-
-#- path to apps
-setenv CFDEM_APP_DIR $FOAM_USER_APPBIN
-
-#------------------------------------------------------------------------------
-#- settings for lpp postproc tool
-#------------------------------------------------------------------------------
-
-#- nr of procs for lpp tool
-setenv CFDEM_LPP_NPROCS 4
-
-#- nr of procs for lpp tool
-setenv CFDEM_LPP_CHUNKSIZE 1
-
-#- shortcut to run lpp
-alias lpp 'python -i $CFDEM_LPP_DIR/lpp.py --cpunum $CFDEM_LPP_NPROCS --chunksize $CFDEM_LPP_CHUNKSIZE \!:1'
-
-#------------------------------------------------------------------------------
-#- aliases for easy navigation (no changes necessary)
-#------------------------------------------------------------------------------
-#- shortcut to cfdem path
-alias cfdem 'cd $CFDEM_PROJECT_DIR'
-
-#- shortcut to src path
-alias cfdemSrc 'cd $CFDEM_SRC_DIR'
-
-#- shortcut to tutorial path
-alias cfdemTut 'cd $CFDEM_TUT_DIR'
-
-#- shortcut to solver path
-alias cfdemSol 'cd $CFDEM_SOLVER_DIR'
-
-#- shortcut to utilities path
-alias cfdemUt 'cd $CFDEM_UT_DIR'
-
-#- shortcut to run path
-alias cfdemRun 'cd $CFDEM_PROJECT_USER_DIR/run'
-
-#- shortcut to user solver path
-alias cfdemUsrSol 'cd $CFDEM_PROJECT_USER_DIR/applications/solvers'
-
-#- shortcut to documentation path
-alias cfdemDoc 'cd $CFDEM_DOC_DIR'
-
-#- shortcut to open the doxygen with firefox
-alias cfdemDox 'firefox $CFDEM_DOC_DIR/doxygen/html/index.html'
-
-#- shortcut to LIGGGHTS path
-alias cfdemLIG 'cd $CFDEM_LIGGGHTS_SRC_DIR'
-
-#- shortcut to system test
-alias cfdemSysTest 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cfdemSystemTest.sh'
-
-#- shortcut to pull LIGGGHTS
-alias cfdemPullLIG 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/pullLIGGGHTS.sh'
-
-#- shortcut to pull CFDEMcoupling
-alias cfdemPullCFDEMcoupling 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/pullCFDEMcoupling.sh'
-
-#- shortcut to clean CFDEM
-alias cfdemCleanCFDEM 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/cleanCFDEMcoupling.sh'
-
-#- shortcut to compile LIGGGHTS + sublibraries
-alias cfdemCompLIG 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileLIGGGHTS.sh'
-
-#- shortcut to compile CFDEMcoupling +LIGGGHTS
-alias cfdemCompCFDEMall 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_all.sh'
-
-#- shortcut to compile CFDEMcoupling (src+solvers)
-alias cfdemCompCFDEM 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling.sh'
-
-#- shortcut to compile CFDEMcoupling src
-alias cfdemCompCFDEMsrc 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_src.sh'
-
-#- shortcut to compile CFDEMcoupling solvers
-alias cfdemCompCFDEMsol 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_sol.sh'
-
-#- shortcut to compile CFDEMcoupling utilities
-alias cfdemCompCFDEMuti 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh'
-
-#- shortcut to test basic tutorials
-alias cfdemTestTUT 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/testTutorials.sh'
-
-#- shortcut to visualize the clock model data
-alias vizClock 'python $CFDEM_UT_DIR/vizClock/matPlot.py'
-
-#- shortcut to run liggghts in serial
-alias cfdemLiggghts '$CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME'
-
-#- shortcut to run liggghts in parallel (no fcts in csh...)
-alias cfdemLiggghtsPar 'echo "mpirun -np xx -machinefile mynodes $CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME < in.liggghts_init"'
-
-#- shortcut to run liggghts in parallel
-# unfortunately no functions available in csh
-
-#- shortcut to run lpp
-alias lpp 'python -i $CFDEM_LPP_DIR/lpp.py \!:1'
-
-# check if the directory exists
-if ( -d "$CFDEM_PROJECT_USER_DIR" ) then
-    :
-else
-    if ( -d "$CFDEM_PROJECT_DIR" ) then
-        echo "make new dirs $CFDEM_PROJECT_USER_DIR ? (y/n)"
-	set YN=$<
-        if ( $YN == "y" ) then
-            mkdir -p $CFDEM_PROJECT_USER_DIR
-            cd $CFDEM_PROJECT_USER_DIR
-            mkdir run
-            mkdir -p log/logFilesCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
-            mkdir -p applications/solvers
-        else
-            echo "aborted by user."
-            exit
-        endif
-    else
-        echo "error in CFDEMcoupling's cshrc."
-        exit
-    endif
-endif

src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.C

@@ -78,55 +78,71 @@ fileName IOModel::createLagrangianDir(fileName path) const
 fileName IOModel::buildFilePath(word dirName) const
 {
     // create file structure
-	fileName path("");
+    fileName path("");
     if(parOutput_)
     {
-    	path=fileName(particleCloud_.mesh().time().path()/particleCloud_.mesh().time().timeName()/dirName/"particleCloud");
-    	mkDir(path,0777);
-    } else
+        path = fileName(particleCloud_.mesh().time().path()/particleCloud_.mesh().time().timeName()/dirName/"particleCloud");
+        mkDir(path,0777);
+    }
+    else
     {
-		path=fileName("."/dirName);
-    	mkDir(path,0777);
-    	mkDir(fileName(path/"constant"),0777);
-    	OFstream* stubFile = new OFstream(fileName(path/"particles.foam"));
-    	delete stubFile;
-        }
+        path = fileName("."/dirName);
+        mkDir(path,0777);
+        mkDir(path/"constant",0777);
+        OFstream stubFile(path/"particles.foam");
+    }
     return path;
 }
 
-void IOModel::streamDataToPath(fileName path, double** array,int nPProc,word name,word type,word className,word finaliser) const
+void IOModel::streamDataToPath(fileName path, double** array,int nPProc,word name,word type,word className) const
 {
-    vector vec;
-    OFstream* fileStream = new OFstream(fileName(path/name));
-    *fileStream << "FoamFile\n";
-    *fileStream << "{version 2.0; format ascii;class "<< className << "; location 0;object  "<< name <<";}\n";
-    *fileStream << nPProc <<"\n";
+    OFstream fileStream(path/name);
 
+    fileStream
+        << "FoamFile\n{\n"
+        << "    version     " << fileStream.version() << ";\n"
+        << "    format      " << fileStream.format() << ";\n"
+        << "    class       " << className << ";\n"
+        << "    location    " << 0 << ";\n"
+        << "    object      " << name << ";\n"
+        << "}" << nl;
 
-    if(type!="origProcId")*fileStream << "(\n";
-    else if(type=="origProcId")
+    fileStream << nPProc <<"\n";
+
+    if (type == "origProcId")
     {
-        if(nPProc>0) *fileStream <<"{0}"<< "\n";
-        else *fileStream <<"{}"<< "\n";
+        if (nPProc > 0) fileStream << "{0}" << "\n";
+        else fileStream << "{}" << "\n";
+        return;
     }
 
-    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    fileStream << token::BEGIN_LIST << nl;
+
+    int ** cellIDs = particleCloud_.cellIDs();
+    for (int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
-        if (particleCloud_.cellIDs()[index][0] > -1) // particle Found
+        if (cellIDs[index][0] > -1) // particle Found
         {
-            if (type=="scalar"){
-                *fileStream << array[index][0] << " \n";
-            }else if (type=="position" || type=="vector"){
-                for(int i=0;i<3;i++) vec[i] = array[index][i];
-                *fileStream <<"( "<< vec[0] <<" "<<vec[1]<<" "<<vec[2]<<" ) "<< finaliser << " \n";
-            }else if (type=="label"){
-                *fileStream << index << finaliser << " \n";
+            if (type == "scalar")
+            {
+                fileStream << array[index][0] << " \n";
+            }
+            else if (type == "position")
+            {
+                fileStream << "( "<< array[index][0] << " " << array[index][1] << " " << array[index][2] << " ) " << cellIDs[index][0] << " \n";
+            }
+            else if (type == "vector")
+            {
+                fileStream << "( "<< array[index][0] << " " << array[index][1] << " " << array[index][2] << " ) " << " \n";
+            }
+            else if (type == "label")
+            {
+                fileStream << index << " \n";
             }
         }
     }
 
-    if(type!="origProcId")*fileStream << ")\n";
-    delete fileStream;
+    fileStream << token::END_LIST << nl;
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -143,12 +159,12 @@ IOModel::IOModel
     time_(sm.mesh().time()),
     parOutput_(true)
 {
-	if (
+    if (
             particleCloud_.dataExchangeM().myType()=="oneWayVTK" ||
             dict_.found("serialOutput")
        )
     {
-        parOutput_=false;
+        parOutput_ = false;
         Warning << "IO model is in serial write mode, only data on proc 0 is written" << endl;
     }
 }

src/lagrangian/cfdemParticle/subModels/IOModel/IOModel/IOModel.H

@@ -121,7 +121,7 @@ public:
 
     fileName buildFilePath(word) const;
 
-    void streamDataToPath(fileName,double**,int,word,word type,word className,word finaliser) const;
+    void streamDataToPath(fileName,double**,int,word,word type,word className) const;
 };
 
 

src/lagrangian/cfdemParticle/subModels/IOModel/basicIO/basicIO.C

@@ -92,19 +92,19 @@ int basicIO::dumpDEMdata() const
         {
             Info << "createTimeDir(path_), path="<<path_ << endl;
             Info << "lagPath_=createTimeDir(fileName(lagPath_/lagrangian)), lagPath="<<path_ << endl;
-        	lagPath_=createTimeDir(path_);
-        	lagPath_=createTimeDir(fileName(lagPath_/"lagrangian"));
+            lagPath_ = createTimeDir(path_);
+            lagPath_ = createTimeDir(fileName(lagPath_/"lagrangian"));
         }
         // calc the number of particles on proc
         int count(0);
         for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
             if (particleCloud_.cellIDs()[index][0] > -1) count++;
-        nPProc_=count;
-        
+        nPProc_ = count;
+
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.positions(),nPProc_,"positions","vector","Cloud<passiveParticle>","0");
-        streamDataToPath(lagPath_, particleCloud_.velocities(),nPProc_,"v","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.radii(),nPProc_,"r","scalar","scalarField","");
+        streamDataToPath(lagPath_, particleCloud_.positions(),nPProc_,"positions","position","Cloud<passiveParticle>");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),nPProc_,"v","vector","vectorField");
+        streamDataToPath(lagPath_, particleCloud_.radii(),nPProc_,"r","scalar","scalarField");
     }
     return nPProc_;
 }

src/lagrangian/cfdemParticle/subModels/IOModel/sophIO/sophIO.C

@@ -84,10 +84,10 @@ int sophIO::dumpDEMdata() const
         npProcs=basicIO::dumpDEMdata();
 
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.voidfractions(),npProcs,"voidfractions","scalar","scalarField","");
-        streamDataToPath(lagPath_, particleCloud_.impForces(),npProcs,"impForces","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.expForces(),npProcs,"expForces","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.DEMForces(),npProcs,"DEMForces","vector","vectorField","");
+        streamDataToPath(lagPath_, particleCloud_.voidfractions(),npProcs,"voidfractions","scalar","scalarField");
+        streamDataToPath(lagPath_, particleCloud_.impForces(),npProcs,"impForces","vector","vectorField");
+        streamDataToPath(lagPath_, particleCloud_.expForces(),npProcs,"expForces","vector","vectorField");
+        streamDataToPath(lagPath_, particleCloud_.DEMForces(),npProcs,"DEMForces","vector","vectorField");
     }
     return npProcs;
 }

src/lagrangian/cfdemParticle/subModels/IOModel/trackIO/trackIO.C

@@ -82,8 +82,8 @@ int trackIO::dumpDEMdata() const
         npProcs = sophIO::dumpDEMdata();
 
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origId","label","labelField","");
-        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origProcId","origProcId","labelField","");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origId","label","labelField");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origProcId","origProcId","labelField");
     }
 
     return npProcs;

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;
@@ -135,7 +135,7 @@ void averagingModel::setVectorSum
 (
     volVectorField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     double**const& mask
 ) const
 {
@@ -169,7 +169,7 @@ void averagingModel::setVectorSumSimple
 (
     volVectorField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     int nP
 ) const
 {
@@ -232,18 +232,16 @@ void averagingModel::setScalarSum
 void averagingModel::setDSauter
 (
     volScalarField& dSauter,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     label myParticleType
 ) const
 {
     label cellI;
-    scalar valueScal;
     scalar weightP;
     scalar radius(-1);
     scalar radiusPow2(-1);
     scalar radiusPow3(-1);
-    scalar volume(-1);
 
     scalar scale_ = particleCloud_.cg(); //scaling of parcel vs. primary particle diameter
     dSauter = 0.0 * dSauter; //set to zero, because we will use it to calc sum(wi*ri^3)
@@ -264,7 +262,7 @@ void averagingModel::setDSauter
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {
         if(myParticleType!=0) //in case a particle type is specified, only consider particles of the right type
-            if(myParticleType != particleCloud_.particleType(index)) continue; 
+            if(myParticleType != particleCloud_.particleType(index)) continue;
 
         radius         = particleCloud_.radii()[index][0] / scale_; //the primary particle diameter
         radiusPow2 = radius*radius;
@@ -305,17 +303,17 @@ 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;
 }
 
 
-void Foam::averagingModel::undoWeightFields(double**const& mask) const
+void averagingModel::undoWeightFields(double**const& mask) const
 {
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {
@@ -328,41 +326,14 @@ void Foam::averagingModel::undoWeightFields(double**const& mask) const
     }
 }
 
-tmp<volVectorField> Foam::averagingModel::UsInterp() const
+tmp<volVectorField> averagingModel::UsInterp() const
 {
-    tmp<volVectorField> tsource
+    const scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
+
+    return tmp<volVectorField>
     (
-        new volVectorField
-        (
-            IOobject
-            (
-                "Us_averagingModel",
-                particleCloud_.mesh().time().timeName(),
-                particleCloud_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            particleCloud_.mesh(),
-            dimensionedVector
-            (
-                "zero",
-                dimensionSet(0, 1, -1, 0, 0),
-                vector::zero
-            )
-        )
+        new volVectorField("Us_averagingModel", (1. - tsf) * UsPrev_ + tsf * UsNext_)
     );
-
-    if (particleCloud_.dataExchangeM().couplingStep() > 1)
-    {
-        tsource() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
-                    + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_;
-    }
-    else
-    {
-        tsource() = UsNext_;
-    }
-
-    return tsource;
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.H

@@ -118,7 +118,7 @@ public:
         (
             volScalarField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const = 0;
@@ -127,7 +127,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const = 0;
@@ -156,7 +156,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             double**const& mask
         ) const;
 
@@ -164,7 +164,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             int np
         ) const;
 
@@ -179,7 +179,7 @@ public:
         void setDSauter
         (
             volScalarField& dSauter,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             label myParticleType = 0  // can be evaluated for a special type
         ) const;

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

@@ -35,7 +35,7 @@ Description
 #include "addToRunTimeSelectionTable.H"
 #include "voidFractionModel.H"
 
-//#include "mpi.h"
+//#include <mpi.h>
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -78,7 +78,7 @@ void dense::setScalarAverage
 (
     volScalarField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     double**const& mask
 ) const
@@ -101,7 +101,7 @@ void dense::setScalarAverage
                 if (cellI >= 0)
                 {
                     valueScal = value[index][0];
-                    weightP = weight[index][0];
+                    weightP = weight[index][subCell];
 
                     // first entry in this cell
                     if(weightField[cellI] == 0)
@@ -127,7 +127,7 @@ void dense::setVectorAverage
 (
     volVectorField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     double**const& mask
 ) const

src/lagrangian/cfdemParticle/subModels/averagingModel/dense/dense.H

@@ -82,7 +82,7 @@ public:
         (
             volScalarField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const;
@@ -91,7 +91,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const;

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

@@ -77,7 +77,7 @@ void dilute::setScalarAverage
 (
     volScalarField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     double**const& mask
 ) const
@@ -115,7 +115,7 @@ void dilute::setVectorAverage
 (
     volVectorField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     double**const& mask
 ) const

src/lagrangian/cfdemParticle/subModels/averagingModel/dilute/dilute.H

@@ -84,7 +84,7 @@ public:
         (
             volScalarField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const;
@@ -93,7 +93,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const;

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

@@ -0,0 +1,61 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "chemistryModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(chemistryModel, 0);
+
+defineRunTimeSelectionTable(chemistryModel, dictionary);
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+chemistryModel::chemistryModel
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    dict_(dict),
+    particleCloud_(sm)
+{}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+chemistryModel::~chemistryModel()
+{}
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/chemistryModel/chemistryModel/chemistryModel.H

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef chemistryModel_H
+#define chemistryModel_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "probeModel.H"
+#include "interpolationCellPoint.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class chemistryModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class chemistryModel
+{
+
+protected:
+
+    // Protected data
+        const dictionary& dict_;
+
+        cfdemCloudEnergy& particleCloud_;
+	
+
+       
+public:
+
+    //- Runtime type information
+    TypeName("chemistryModel");
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            chemistryModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                cfdemCloudEnergy& sm
+            ),
+            (dict,sm)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        chemistryModel
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm            
+        );
+
+
+    // Destructor
+
+        virtual ~chemistryModel();
+
+
+    // Selector
+
+        static autoPtr<chemistryModel> New
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Member Functions
+
+       virtual void execute() = 0;
+       
+//        virtual tmp<fvScalarMatrix> Smi (const label i) const = 0;
+//        
+//        virtual tmp<fvScalarMatrix> Sm () const = 0;
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "chemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+autoPtr<chemistryModel> chemistryModel::New
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+{
+    word chemistryModelType
+    (
+        dict.lookup("chemistryModel")
+    );
+    Info<< "Selecting chemistryModel "
+         << chemistryModelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(chemistryModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalError
+            << "chemistryModel::New(const dictionary&, cfdemCloudEnergy&) : "
+            << endl
+            << "    unknown chemistryModelType type "
+            << chemistryModelType
+            << ", constructor not in hash table" << endl << endl
+            << "    Valid chemistryModel types are :"
+            << endl;
+        Info<< dictionaryConstructorTablePtr_->toc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<chemistryModel>(cstrIter()(dict,sm));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "noChemistry.H"
+#include "addToRunTimeSelectionTable.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(noChemistry, 0);
+
+addToRunTimeSelectionTable(chemistryModel, noChemistry, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+noChemistry::noChemistry
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    chemistryModel(dict,sm)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+noChemistry::~noChemistry()
+{}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void noChemistry::execute()
+{}
+
+//tmp<Foam::fvScalarMatrix> noChemistry::Smi(const label i) const
+//{
+//    return tmp<fvScalarMatrix>(new fvScalarMatrix(0, dimMass/dimTime)); 
+//}
+
+//tmp<Foam::fvScalarMatrix> noChemistry::Sm() const
+//{
+//    return tmp<fvScalarMatrix>(new fvScalarMatrix(0, dimMass/dimTime)); 
+//}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/chemistryModel/noChemistry/noChemistry.H

@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+    
+    Description
+    Chemistry turned off
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noChemistry_H
+#define noChemistry_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "chemistryModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class noChemistry Declaration
+\*---------------------------------------------------------------------------*/
+
+
+class noChemistry
+:
+    public chemistryModel
+{
+private:
+
+       
+public:
+
+    //- Runtime type information
+    TypeName("off");
+
+    // Constructors
+
+        //- Construct from components
+        noChemistry
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~noChemistry();
+
+
+    // Member Functions
+
+    void execute();
+    
+//     tmp<fvScalarMatrix> Smi(const label i) const;
+//     
+//     tmp<fvScalarMatrix> Sm() const;
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+                        M.Efe Kinaci, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "species.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "dataExchangeModel.H"
+#include "IFstream.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(species, 0);
+
+addToRunTimeSelectionTable
+(
+        chemistryModel,
+        species,
+        dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+species::species
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    chemistryModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
+    mesh_(sm.mesh()),
+    // define a file name in the coupling properties that contains the species
+    specDict_
+    (
+        IFstream
+        (
+            fileName(propsDict_.lookup("ChemistryFile")).expand()
+        )()
+    ),
+    // create a list from the Species table in the specified species dictionary
+    speciesNames_(specDict_.lookup("species")),
+    mod_spec_names_(speciesNames_.size()),
+    Y_(speciesNames_.size()),                           //volumeScalarFields created in the ts folders
+    concentrations_(speciesNames_.size(),NULL),         //the value of species concentration for every species
+    changeOfSpeciesMass_(speciesNames_.size(),NULL),    //the values that are received from DEM with the name of Modified_+species name
+    changeOfSpeciesMassFields_(speciesNames_.size()),   //the scalar fields generated with the values from Modified_+species names
+    changeOfGasMassField_                               //the total change of Gas Mass field (when the Modified species
+
+    (
+        IOobject
+        (
+            "changeOfGasMassField",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+        mesh_,
+        dimensionedScalar("zero",dimMass/(dimVol*dimTime),0.0)
+    ),
+    tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
+    tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
+    partTempName_(propsDict_.lookupOrDefault<word>("partTempName","partTemp")),
+    partTemp_(NULL),
+    densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
+    rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
+    partRhoName_(propsDict_.lookupOrDefault<word>("partRhoName","partRho")),
+    partRho_(NULL)
+  // voidfraction and velocity fields can be included by wish
+  /*  voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
+      voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+      velFieldName_(propsDict_.lookup("velFieldName")),
+      U_(sm.mesh().lookup<volVectorField> (velFieldName_)),*/
+
+{
+    Info << " Read species list from: " << specDict_.name() << endl;
+    Info << " Reading species list: " << speciesNames_ << endl;
+
+    for (int i=0; i<speciesNames_.size(); i++)
+    {
+        // Defining the Species volume scalar fields
+        Info << " Looking up species fields " << speciesNames_[i] << endl;
+        volScalarField& Y = const_cast<volScalarField&>
+                (sm.mesh().lookupObject<volScalarField>(speciesNames_[i]));
+        Y_.set(i, &Y);
+
+        // define the modified species names
+        mod_spec_names_[i] = "Modified_" + speciesNames_[i];
+
+        // Check if mod species are correct
+        Info << "Modified species names are: " << mod_spec_names_[i] << endl;
+
+        // Create new volScalarFields for the changed values of the species mass fields
+        changeOfSpeciesMassFields_.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                "ModSpeciesMassField_"+Y_[i].name(),
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("0",dimMass/(dimVol*dimTime), 0)
+            )
+         );
+
+        Info << "what are the concentration names (Y_i): " << Y_[i].name() << endl;
+
+    }
+    allocateMyArrays();
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+species::~species()
+{
+    particleCloud_.dataExchangeM().destroy(partTemp_,1);
+    particleCloud_.dataExchangeM().destroy(partRho_,1);
+
+    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(concentrations_[i],1);
+    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(changeOfSpeciesMass_[i],1);
+
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+void species::allocateMyArrays() const
+{
+  // could be implemented similarly as forcemodel LaEuScalarTemp
+
+    // get memory for 2d arrays
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partRho_,initVal,1);
+    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);
+
+    for (int i=0; i<speciesNames_.size(); i++)
+    {
+        particleCloud_.dataExchangeM().allocateArray(concentrations_[i],initVal,1);
+        particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1);
+    }
+}
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void species::execute()
+{
+    // realloc the arrays
+    allocateMyArrays();
+
+    // get Y_i, T, rho at particle positions, fill arrays with them and push to LIGGGHTS
+
+    label  cellI=0;
+    scalar Tfluid(0);
+    scalar rhofluid(0);
+    List<scalar> Yfluid_;
+    Yfluid_.setSize(speciesNames_.size());
+
+    // defining interpolators for T, rho
+    interpolationCellPoint <scalar> TInterpolator_(tempField_);
+    interpolationCellPoint <scalar> rhoInterpolator_(rho_);
+
+    for (int index=0; index<particleCloud_.numberOfParticles(); index++)
+    {
+        cellI=particleCloud_.cellIDs()[index][0];
+        if (cellI >= 0)
+        {
+            if(interpolation_)
+            {
+                vector position     =   particleCloud_.position(index);
+                Tfluid              =   TInterpolator_.interpolate(position,cellI);
+                rhofluid            =   rhoInterpolator_.interpolate(position,cellI);
+            }
+            else
+            {
+                Tfluid = tempField_[cellI];
+                rhofluid=rho_[cellI];
+            }
+
+            //fill arrays
+            partTemp_[index][0]=Tfluid;
+            partRho_[index][0]=rhofluid;
+            for (int i=0; i<speciesNames_.size();i++)
+            {
+                Yfluid_[i] = Y_[i][cellI];
+                concentrations_[i][index][0] = Yfluid_[i];
+            }
+        }
+
+        if(particleCloud_.verbose() && index >= 0 && index < 2)
+        {
+            /*for(int i =0; i<speciesNames_.size();i++)
+            {
+                Info << "Y_i = " << Y_[i].name() << endl;
+                Info << "concentrations = " << concentrations_[i][index][0] << endl;
+                Info << "partRho_[index][0] = " << partRho_[index][0] << endl;
+                Info << "rhofluid =" << rhofluid << endl;
+                Info << "Yfluid = " << Yfluid_[i] << endl;
+                Info << "partTemp_[index][0] = " << partTemp_[index][0] << endl;
+                Info << "Tfluid = " << Tfluid << endl  ;
+            }*/
+        }
+    }
+
+        // give DEM data
+        particleCloud_.dataExchangeM().giveData(partTempName_, "scalar-atom", partTemp_);
+        particleCloud_.dataExchangeM().giveData(partRhoName_, "scalar-atom", partRho_);
+        for (int i=0; i<speciesNames_.size();i++)
+        {
+            particleCloud_.dataExchangeM().giveData(Y_[i].name(),"scalar-atom",concentrations_[i]);
+        };
+
+        Info << "give data done" << endl;
+
+
+
+        // pull changeOfSpeciesMass_, transform onto fields changeOfSpeciesMassFields_, add them up on changeOfGasMassField_
+        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].primitiveFieldRef() = 0.0;
+            changeOfSpeciesMassFields_[i].boundaryFieldRef() = 0.0;
+            particleCloud_.averagingM().setScalarSum
+            (
+                changeOfSpeciesMassFields_[i],
+                changeOfSpeciesMass_[i],
+                particleCloud_.particleWeights(),
+                NULL
+            );
+
+            // take care for implementation in LIGGGHTS: species produced from particles defined positive
+            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;
+        }
+        Info << "get data done" << endl;
+
+}
+
+/*tmp<Foam::fvScalarMatrix> species::Smi(const label i) const
+{
+    return tmp<fvScalarMatrix>(new fvScalarMatrix(changeOfSpeciesMassFields_[i], dimMass/dimTime));
+}
+
+tmp<Foam::fvScalarMatrix> species::Sm() const
+{
+    return tmp<fvScalarMatrix>(new fvScalarMatrix(changeOfGasMassField_, dimMass/dimTime));
+}*/
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+    
+    Description
+    Model to communicate species concentration and changes due to reactions
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef species_H
+#define species_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "chemistryModel.H"
+
+#include "HashPtrTable.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class species Declaration
+\*---------------------------------------------------------------------------*/
+
+// for future use:
+// + communicate every N steps
+
+class species
+:
+    public chemistryModel
+{
+private:
+
+    dictionary propsDict_;
+
+    bool interpolation_;
+
+    const fvMesh& mesh_;
+
+    // Species dictionary - To be read where the species are from
+    dictionary specDict_;
+
+    wordList speciesNames_;
+    
+    wordList mod_spec_names_;
+
+    PtrList<volScalarField> Y_;
+
+    mutable List<double**> concentrations_;
+    
+    mutable List<double**> changeOfSpeciesMass_;
+
+    PtrList<volScalarField> changeOfSpeciesMassFields_;
+    
+    volScalarField changeOfGasMassField_;
+    
+    word tempFieldName_;
+
+    const volScalarField& tempField_;              // ref to gas temperature field
+
+    word partTempName_;
+
+    mutable double **partTemp_;                    // gas temperature at particle positions
+
+    word densityFieldName_;
+
+    const volScalarField& rho_;
+
+    word partRhoName_;
+    
+    mutable double **partRho_;                     // gas density at particle positions
+    
+    //word velFieldName_;                          // when velocity is to be calculated in the species source code
+
+    //const volVectorField& U_;
+
+    //word voidfractionFieldName_;                 // when voidfraction is to be calculated in the species source code
+
+    //const volScalarField& voidfraction_;
+
+
+    void allocateMyArrays() const;
+       
+public:
+
+    //- Runtime type information
+    TypeName("species");
+
+    // Constructors
+
+        //- Construct from components
+        species
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~species();
+
+
+    // Member Functions
+
+    void execute();
+
+    /*tmp<fvScalarMatrix> Smi(const label i) const;
+
+    tmp<fvScalarMatrix> Sm() const;*/
+    
+    
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/clockModel.C

@@ -30,10 +30,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
-#include "mpi.h"
+#include <mpi.h>
 #include "clockModel.H"
 #include <unistd.h>
-
+#include <time.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,475 +48,382 @@ defineRunTimeSelectionTable(clockModel, dictionary);
 
 // * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * //
 
-void Foam::clockModel::start(int pos) const
+void clockModel::start(int pos) const
 {
-	start(pos,"");
-	return;
+    start(pos,"");
 }
 
-void Foam::clockModel::start(int pos,std::string ident) const
+void clockModel::start(int pos, const std::string& ident) const
 {
-    if(particleCloud_.mesh().time().value() > startTime_)
+    if (time_.value() > startTime_)
     {
-    	if (pos >= n_) // alternatively one fixed size?
-    	{
-    		n_ = 2*n_;
-	    	deltaT_.resize(n_,0);
-	    	identifier_.resize(n_,"");
-	    	nOfRuns_.resize(n_,0);
-	    	level_.resize(n_,-1);
-	    	parent_.resize(n_,-2);
-	    }
-	    identifier_[pos]=ident;
-	    level_[pos] = curLev_;
-	    curLev_ += 1;
-	    parent_[pos]=curParent_;
-	    curParent_ = pos;
-	    nOfRuns_[pos] += 1;
-	    deltaT_[pos]-=std::clock();
+        if (pos >= n_) // alternatively one fixed size?
+        {
+            n_ += 32;
+            deltaT_.resize(n_,0);
+            identifier_.resize(n_,"");
+            nOfRuns_.resize(n_,0);
+            level_.resize(n_,-1);
+            parent_.resize(n_,-2);
+        }
+        identifier_[pos]=ident;
+        level_[pos] = curLev_;
+        curLev_ += 1;
+        parent_[pos]=curParent_;
+        curParent_ = pos;
+        nOfRuns_[pos] += 1;
+        deltaT_[pos] -= std::clock();
     }
-	return;
 }
 
-void Foam::clockModel::stop() const
+void clockModel::stop() const
 {
-    if(particleCloud_.mesh().time().value() > startTime_)
+    if (time_.value() > startTime_)
     {
-    	deltaT_[curParent_]+=std::clock();
-	    curLev_ -= 1;
-	    if (curParent_ >= 0)
-	    {
-	    	curParent_ = parent_[curParent_];
-	    }
-	    else
-	    {
-	    	curParent_ = -1;
-	    }
+        deltaT_[curParent_] += std::clock();
+        curLev_ -= 1;
+        if (curParent_ >= 0)
+        {
+            curParent_ = parent_[curParent_];
+        }
+        else
+        {
+            curParent_ = -1;
+        }
     }
-	return;
 }
 
-void Foam::clockModel::stop(std::string ident) const
+void clockModel::stop(const std::string& ident) const
 {
-    if(particleCloud_.mesh().time().value() > startTime_)
+    if (time_.value() > startTime_)
     {
-    	deltaT_[curParent_] += std::clock();
-	    if (curParent_ > 0 && identifier_[curParent_].compare(ident)!=0)
-	    {
-	    	Pout<<"Warning: stop identifier did not equal start identifier! "<<ident<<" & "<<identifier_[curParent_]<<nl;
-	    }
-	    curLev_ -= 1;
-	    if (curParent_ >= 0)
-	    {
-	    	curParent_ = parent_[curParent_];
-	    }
-	    else
-	    {
-	    	curParent_ = -1;
-	    }
+        deltaT_[curParent_] += std::clock();
+        if (curParent_ > 0 && identifier_[curParent_].compare(ident) != 0)
+        {
+            Pout << "Warning: stop identifier did not equal start identifier! " << ident << " & " << identifier_[curParent_] << nl;
+        }
+        curLev_ -= 1;
+        if (curParent_ >= 0)
+        {
+            curParent_ = parent_[curParent_];
+        }
+        else
+        {
+            curParent_ = -1;
+        }
     }
-	return;
 }
 
-std::string Foam::clockModel::eval() const
+std::string clockModel::eval() const
 {
-	std::string msg = "Measurements in CPU-seconds:";
-	msg.append("\n");
-	msg.append("Name \t deltaT \t nOfRuns \t level \t parentNr \t parentName \n");
-	std::ostringstream strs;
-	strs.setf(std::ios_base::scientific);
-	std::vector<int> shifts = calcShift();
+    std::ostringstream strs("Measurements in CPU-seconds:\n");
+    strs << "Name\tdeltaT\tnOfRuns\tlevel\tparentNr\tparentName\n";
+    strs.setf(std::ios_base::scientific);
+    std::vector<int> shifts = calcShift();
 
-	for (int i=0;i<n_;i++)
-	{
-		if (parent_[i] != -2)
-		{
-			msg.append(identifier_[i]);
-			msg.append("\t");
-			strs << (double(deltaT_[i])/(CLOCKS_PER_SEC));
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
+    for (int i=0; i<n_; i++)
+    {
+        if (parent_[i] != -2)
+        {
+            strs << identifier_[i] << "\t";
+            strs << static_cast<double>(deltaT_[i])/(CLOCKS_PER_SEC) << "\t";
+            strs << nOfRuns_[i] << "\t";
+            strs << level_[i] << "\t";
 
-			strs << nOfRuns_[i];
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
-
-			strs << level_[i];
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
-
-			if (parent_[i] >= 0)
-			{
-				strs << (shifts[parent_[i]]);
-			}
-			else
-			{
-				strs << parent_[i];
-			}
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
-
-			if (parent_[i] >= 0)
-			{
-				strs << identifier_[parent_[i]];
-			}
-			else
-			{
-				strs << "none";
-			}
-
-			msg.append(strs.str());
-			msg.append("\n");
-			strs.str("");
-		}
-	}
-	return msg;
+            if (parent_[i] >= 0)
+            {
+                strs << (shifts[parent_[i]]) << "\t";
+                strs << identifier_[parent_[i]] << "\n";
+            }
+            else
+            {
+                strs << parent_[i] << "\t";
+                strs << "none\n";
+            }
+        }
+    }
+    return strs.str();
 }
 
-void Foam::clockModel::evalFile() const
+void clockModel::evalFile() const
 {
-	std::ofstream outFile;
-	std::string fileName(path_/"timeEval.txt");
-	outFile.open(fileName.data(),ios_base::trunc);
-	outFile << "Time Evaluation"<<nl;
-	outFile << eval();
-	outFile.close();
-
+    std::ofstream outFile;
+    std::string fileName(path_/"timeEval.txt");
+    outFile.open(fileName.c_str(), ios_base::trunc);
+    outFile << "Time Evaluation" << nl;
+    outFile << eval();
+    outFile.close();
 }
 
-void Foam::clockModel::evalPar() const
+void clockModel::evalPar() const
 {
-	int myrank=-10;
-	MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
-	int numprocs=-10;
-	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
-	std::ofstream outFile;
-	std::ostringstream strs;
-	strs.setf(std::ios_base::scientific);
-	std::string fileName(path_/"timeEval_");
-	strs << myrank;
-	fileName.append(strs.str());
-	fileName.append(".txt");
-	outFile.open(fileName.data(),ios_base::trunc);
-	outFile << "Time Evaluation for Processor Nr." << myrank <<nl;
-	outFile << eval();
-	outFile.close();
+    int myrank, numprocs;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
 
-	// MPI_REDUCE SUM NODES
-	MPI_Barrier(MPI_COMM_WORLD);
-	strs.str("");
-	std::string msg = "Parallel Measurements in CPU-seconds of all Processors (starting after first t.s.):";
-	msg.append("\n");
-	msg.append("Name \t avgdeltaT \t maxdeltaT \t nOfRuns \t level \t parentNr \t parentName \n");
-	double buffOut=0.;
-    double buffIn=0.;
-	std::vector<int> shifts = calcShift();
+    std::ofstream outFile;
+    std::ostringstream strs;
+    strs.setf(std::ios_base::scientific);
 
-	for (int i=0;i<n_;i++)
-	{
-		if (parent_[i]!=-2)
-		{
-			msg.append(identifier_[i]);
-			msg.append("\t");
+    std::string fileName(path_/"timeEval_");
+    strs << myrank << ".txt";
+    fileName.append(strs.str());
 
-			buffIn = double(deltaT_[i])/(CLOCKS_PER_SEC);
-			MPI_Reduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM,0 , MPI_COMM_WORLD);
-			MPI_Barrier(MPI_COMM_WORLD);
-			strs << (buffOut)/(numprocs);
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
+    outFile.open(fileName.c_str(), ios_base::trunc);
+    outFile << "Time Evaluation for Processor Nr." << myrank << nl;
+    outFile << eval();
+    outFile.close();
 
-			//buffIn = double(deltaT[i])/(CLOCKS_PER_SEC);
-			MPI_Reduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_MAX,0 , MPI_COMM_WORLD);
-			MPI_Barrier(MPI_COMM_WORLD);
-			strs << (buffOut);
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
+    // MPI_REDUCE SUM NODES
+    MPI_Barrier(MPI_COMM_WORLD);
+    strs.str("Parallel Measurements in CPU-seconds of all Processors (starting after first t.s.):\n");
+    strs << "Name\tavgdeltaT\tmaxdeltaT\tnOfRuns\tlevel\tparentNr\tparentName\n";
+    double buffOut = 0.;
+    double buffIn = 0.;
+    std::vector<int> shifts = calcShift();
 
-			strs << nOfRuns_[i];
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
+    for (int i=0; i<n_; i++)
+    {
+        if (parent_[i] != -2)
+        {
+            strs << identifier_[i] << "\t";
 
-			strs << level_[i];
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
+            buffIn = static_cast<double>(deltaT_[i])/(CLOCKS_PER_SEC);
 
-			if (parent_[i] >= 0)
-			{
-				strs << (shifts[parent_[i]]);
-			}
-			else
-			{
-				strs << parent_[i];
-			}
-			msg.append(strs.str());
-			msg.append("\t");
-			strs.str("");
+            MPI_Reduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+            MPI_Barrier(MPI_COMM_WORLD);
 
-			if (parent_[i] >= 0)
-			{
-				strs << identifier_[parent_[i]];
-			}
-			else
-			{
-				strs << "none";
-			}
+            strs << buffOut/numprocs << "\t";
 
-			msg.append(strs.str());
-			msg.append("\n");
-			strs.str("");
-		}
-	}
-	MPI_Barrier(MPI_COMM_WORLD);
+            MPI_Reduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
+            MPI_Barrier(MPI_COMM_WORLD);
 
-	if(myrank == 0)
-	{
+            strs << buffOut << "\t";
+            strs << nOfRuns_[i] << "\t";
+            strs << level_[i] << "\t";
+
+            if (parent_[i] >= 0)
+            {
+                strs << (shifts[parent_[i]]) << "\t";
+                strs << identifier_[parent_[i]] << "\n";
+            }
+            else
+            {
+                strs << parent_[i] << "\t";
+                strs << "none\n";
+            }
+        }
+    }
+    MPI_Barrier(MPI_COMM_WORLD);
+
+    if (myrank == 0)
+    {
         std::string fileName(path_/"timeEvalFull.txt");
-		outFile.open(fileName.data(),ios_base::trunc);
-		outFile << msg;
-		outFile.close();
-	}
-
-	return;
+        outFile.open(fileName.data(),ios_base::trunc);
+        outFile << strs.str();
+        outFile.close();
+    }
 }
 
 
-void Foam::clockModel::initElems()
+void clockModel::initElems()
 {
     //init elems
-    for (int i = 0;i < n_; i++)
+    for (int i = 0; i < n_; ++i)
     {
-    	deltaT_[i] = 0;
-    	identifier_[i] = "";
-    	nOfRuns_[i] = 0;
-    	level_[i] = -1;
-    	parent_[i] = -2;
+        deltaT_[i] = 0;
+        identifier_[i].clear();
+        nOfRuns_[i] = 0;
+        level_[i] = -1;
+        parent_[i] = -2;
     }
 }
 
-std::vector<int> Foam::clockModel::calcShift() const
+std::vector<int> clockModel::calcShift() const
 {
-	std::vector<int> shifts = std::vector<int> (n_);
-	shifts[0]=0;
-	for (int i=1;i<n_;i++)
-	{
-		if (parent_[i] == -2)
-		{
-			shifts[i] = shifts[i-1];
-		}
-		else
-		{
-			shifts[i] = shifts[i-1]+1;
-		}
-	}
-	return shifts;
+    std::vector<int> shifts(n_, 0);
+
+    for (int i=1; i<n_; ++i)
+    {
+        if (parent_[i] == -2)
+        {
+            shifts[i] = shifts[i-1];
+        }
+        else
+        {
+            shifts[i] = shifts[i-1]+1;
+        }
+    }
+    return shifts;
 }
 
-void Foam::clockModel::normHist() const
+void clockModel::normHist() const
 {
-	int myrank=-10;
-	MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
-	int numprocs=-10;
-	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
-	double buffOut=0.;
+    int myrank, numprocs;
+    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
+    double buffOut = 0.;
     double buffIn=0.;
-	
-	Info << "==========================" << endl;
+
+    Info << "==========================" << endl;
     Info << " PROCESSOR LOAD HISTOGRAM" << endl;
 
-	//Global = 1
-	buffIn = double(deltaT_[1]);
-	MPI_Allreduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-    if(buffOut>SMALL) buffIn /= buffOut;
-    plotHist(buffIn,identifier_[1],numprocs,myrank);
+    //Global = 1
+    buffIn = double(deltaT_[1]);
+    MPI_Allreduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    if (buffOut > SMALL) buffIn /= buffOut;
+    plotHist(buffIn, identifier_[1], numprocs, myrank);
 
-	//LIGGGHTS = 3
-	buffIn = double(deltaT_[3]);
-	MPI_Allreduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-    if(buffOut>SMALL) buffIn /= buffOut;
-    plotHist(buffIn,identifier_[3],numprocs,myrank);
+    //LIGGGHTS = 3
+    buffIn = double(deltaT_[3]);
+    MPI_Allreduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    if (buffOut > SMALL) buffIn /= buffOut;
+    plotHist(buffIn, identifier_[3], numprocs, myrank);
 
-	//Coupling - LIGGGHTS = 2 - 3
-	buffIn = double(deltaT_[2]) - buffIn;
-	MPI_Allreduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-    if(buffOut>SMALL) buffIn /= buffOut;
-    plotHist(buffIn,"Coupling (routines)",numprocs,myrank);
+    //Coupling - LIGGGHTS = 2 - 3
+    buffIn = double(deltaT_[2]) - buffIn;
+    MPI_Allreduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    if (buffOut > SMALL) buffIn /= buffOut;
+    plotHist(buffIn, "Coupling (routines)", numprocs, myrank);
 
-	//Flow = 26
-	buffIn = double(deltaT_[26]);
-	MPI_Allreduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-    if(buffOut>SMALL) buffIn /= buffOut;
-    plotHist(buffIn,identifier_[26],numprocs,myrank);
+    //Flow = 26
+    buffIn = double(deltaT_[26]);
+    MPI_Allreduce(&buffIn, &buffOut, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    if (buffOut > SMALL) buffIn /= buffOut;
+    plotHist(buffIn, identifier_[26], numprocs, myrank);
     Info << "===========================" << endl;
 
-	getRAMUsage();
-	return;
+    getRAMUsage();
 }
 
-void Foam::clockModel::plotHist(double buffIn,std::string identifier,int numprocs,int myrank) const
+void clockModel::plotHist(double buffIn,const std::string& identifier,int numprocs,int myrank) const
 {
-/*  // version using double*, problem: no alloc for double * and MPI
-    double* globalTime=NULL;
-    double* globalTime_all=NULL;
-    particleCloud_.dataExchangeM().allocateArray(globalTime,0.,numprocs);
-    particleCloud_.dataExchangeM().allocateArray(globalTime_all,0.,numprocs);  
-   
-    globalTime[myrank]=buffIn;
-    MPI_Allreduce(globalTime, globalTime_all, numprocs, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    double* globalTime_all = NULL;
 
-    if(myrank==0)
-        for(int j=0;j<numprocs;j++)
-            printf("%4f  ",globalTime_all[j]);
-    Info << "\t" <<identifier << endl;
+    if (myrank == 0) globalTime_all = new double[numprocs];
+    MPI_Gather(&buffIn, 1, MPI_DOUBLE, globalTime_all, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
 
-    particleCloud_.dataExchangeM().destroy(globalTime);
-    particleCloud_.dataExchangeM().destroy(globalTime_all);*/
+    if (myrank == 0)
+        for (int j=0; j<numprocs; j++)
+            printf("%4f  ", globalTime_all[j]);
 
+    Info << "\t" << identifier << endl;
 
-    double** globalTime=NULL;
-    double** globalTime_all=NULL;
-    particleCloud_.dataExchangeM().allocateArray(globalTime,0.,1,numprocs);
-    particleCloud_.dataExchangeM().allocateArray(globalTime_all,0.,1,numprocs);  
-   
-    globalTime[0][myrank]=buffIn;
-    MPI_Allreduce(globalTime[0], globalTime_all[0], numprocs, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-
-    if(myrank==0)
-        for(int j=0;j<numprocs;j++)
-            printf("%4f  ",globalTime_all[0][j]);
-    Info << "\t" <<identifier << endl;
-
-    particleCloud_.dataExchangeM().destroy(globalTime,1);
-    particleCloud_.dataExchangeM().destroy(globalTime_all,1);
+    delete [] globalTime_all;
 }
 
-void Foam::clockModel::Hist() const
+void clockModel::Hist() const
 {
-    int myrank=-10;
-    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+    int myrank = -1;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
 
     //Global = 1 / Coupling = 2 / LIGGGHTS = 3 /Flow = 26
 
     //Global = 1
-	Pout << "[" << myrank << "]: " << identifier_[1] << " " << (deltaT_[1]/CLOCKS_PER_SEC) << '\n';
-	//LIGGGHTS = 3
-	Pout << "[" << myrank << "]: " << identifier_[3] << " " << (deltaT_[3]/CLOCKS_PER_SEC) << '\n';
-	//Coupling - LIGGGHTS = 2 - 3
-	Pout << "[" << myrank << "]: " << "Coupling - LIGGGHTS" << " " << ((deltaT_[2]-deltaT_[3])/CLOCKS_PER_SEC) << '\n';
-	//Flow = 26
-	Pout << "[" << myrank << "]: " << identifier_[26] << " " << (deltaT_[26]/CLOCKS_PER_SEC) << '\n';
-		
-	return;
+    Pout << "[" << myrank << "]: " << identifier_[1] << " " << (deltaT_[1]/CLOCKS_PER_SEC) << '\n';
+    //LIGGGHTS = 3
+    Pout << "[" << myrank << "]: " << identifier_[3] << " " << (deltaT_[3]/CLOCKS_PER_SEC) << '\n';
+    //Coupling - LIGGGHTS = 2 - 3
+    Pout << "[" << myrank << "]: " << "Coupling - LIGGGHTS" << " " << ((deltaT_[2]-deltaT_[3])/CLOCKS_PER_SEC) << '\n';
+    //Flow = 26
+    Pout << "[" << myrank << "]: " << identifier_[26] << " " << (deltaT_[26]/CLOCKS_PER_SEC) << '\n';
 }
 
-void Foam::clockModel::getRAMUsage() const
+void clockModel::getRAMUsage() const
 {
-	int myrank=-10;
-    	MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
-	int numprocs=-10;
-	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
-	
-	pid_t myPID = getpid();	//get PID of running process
-	//Pout << myPID << "\n";
-	
-	std::string fileName = "/proc/"; //build path to /proc/PID/smaps and open file
-	std::stringstream strs;
-	strs << myPID;
+    int myrank, numprocs;
+    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
 
-	fileName.append(strs.str());
-	fileName.append("/smaps");
-	std::ifstream inFile;
-	inFile.open(fileName.data(),ios_base::in);
+    pid_t myPID = getpid();	//get PID of running process
+    //Pout << myPID << "\n";
 
-	std::string line;
-	int RssMem = 0;
-	int SwapMem = 0;
-	int temp = 0;
-	strs.str("");
-	if (inFile.is_open())	//search in File smaps for Rss and Swap entries
-	{
-		while(inFile.good())
-		{
-			getline(inFile,line);
-			strs.str("");
-			if (line.substr(0,4).compare("Rss:") == 0)
-			{
-			strs << line;
-			strs >> line >> temp;
-			RssMem = RssMem + temp;
-			//Pout << temp << " ";
-			}
-			else if (line.substr(0,5).compare("Swap:") == 0)
-			{
-			strs << line;
-			strs >> line >> temp;
-			SwapMem = SwapMem + temp;
-			//Pout << strs.str() << " ";
-			}
-		
-		}
-	}
-	double SwapMB = (double)SwapMem/1024.0; //kB -> MB
-	double RssMB = (double)RssMem/1024.0;
+    std::string fileName = "/proc/"; //build path to /proc/PID/smaps and open file
+    std::stringstream strs;
+    strs << myPID;
 
-	inFile.close();
+    fileName.append(strs.str());
+    fileName.append("/smaps");
+    std::ifstream inFile;
+    inFile.open(fileName.data(),ios_base::in);
 
-	// set up communication between Procs and plot Stuff
-	Info << " RAM USAGE HISTOGRAM in MB" << endl;
-	plotHist(RssMB,"RSS memory used",numprocs,myrank);
-	if (SwapMem > 0) 
-	{
-		plotHist(SwapMB,"WARNING: Swap",numprocs,myrank);
-	}
-	Info << "===========================" << endl;
+    std::string line;
+    int RssMem = 0;
+    int SwapMem = 0;
+    int temp = 0;
+    strs.str("");
 
-	//Pout << "SWAP Memory used: " << SwapMem <<"MB\n";
-	//Pout << "Rss Memory used: " << RssMem <<"MB\n";
-	
+    if (inFile.is_open())	//search in File smaps for Rss and Swap entries
+    {
+        while (inFile.good())
+        {
+            getline(inFile,line);
+            strs.str("");
+            if (line.substr(0,4).compare("Rss:") == 0)
+            {
+                strs << line;
+                strs >> line >> temp;
+                RssMem = RssMem + temp;
+                //Pout << temp << " ";
+            }
+            else if (line.substr(0,5).compare("Swap:") == 0)
+            {
+                strs << line;
+                strs >> line >> temp;
+                SwapMem = SwapMem + temp;
+                //Pout << strs.str() << " ";
+            }
+        }
+    }
+    double SwapMB = SwapMem/1024.0; //kB -> MB
+    double RssMB = RssMem/1024.0;
 
-	return;
+    inFile.close();
+
+    // set up communication between Procs and plot Stuff
+    Info << " RAM USAGE HISTOGRAM in MB" << endl;
+    plotHist(RssMB,"RSS memory used",numprocs,myrank);
+
+    if (SwapMem > 0)
+    {
+        plotHist(SwapMB,"WARNING: Swap",numprocs,myrank);
+    }
+    Info << "===========================" << endl;
+
+    //Pout << "SWAP Memory used: " << SwapMem <<"MB\n";
+    //Pout << "Rss Memory used: " << RssMem <<"MB\n";
 }
 // * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * //
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::clockModel::clockModel
+clockModel::clockModel
 (
     const dictionary& dict,
-    cfdemCloud& sm
+    const Time& time
 )
 :
     dict_(dict),
-    particleCloud_(sm),
+    time_(time),
     path_("clockData"),
-    startTime_(sm.mesh().time().startTime().value()+sm.mesh().time().deltaT().value()+SMALL),  // delay start of measurement by deltaT
-    //startTime_(0),                                //no delay
-    n_(30),
-    deltaT_(std::vector<clock_t> (n_)),
-    identifier_(std::vector<std::string> (n_)),
-    nOfRuns_(std::vector<int> (n_)),
-    level_(std::vector<short> (n_)),
+    startTime_(time.startTime().value()+time.deltaT().value()+SMALL),  // delay start of measurement by deltaT
+    n_(32),
+    deltaT_(n_),
+    identifier_(n_),
+    nOfRuns_(n_),
+    level_(n_),
     curLev_(0),
-    parent_(std::vector<int> (n_)),
+    parent_(n_),
     curParent_(0)
 {
-
-    Info << "start clock measurement at t >"  << startTime_ << endl;
+    Info << "start clock measurement at t >" << startTime_ << endl;
 }
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::clockModel::~clockModel()
+clockModel::~clockModel()
 {}