Merge pull request #15 from ParticulateFlow/release

Release 17.02

.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,31 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude
+
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lfvOptions \
-    -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,22 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -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,23 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I ../cfdemSolverPiso \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I../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,23 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I../cfdemSolverPiso \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -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

@@ -25,7 +25,7 @@ License
     along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cfdemSolverPisoScalar
+    cfdemSolverPiso
 
 Description
     Transient solver for incompressible flow.
@@ -36,10 +36,13 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
+#include "clockModel.H"
 #include "smoothingModel.H"
 #include "forceModel.H"
 
@@ -50,7 +53,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     // create cfdemCloud
@@ -62,12 +67,14 @@ int main(int argc, char *argv[])
     Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
+        particleCloud.clockM().start(1,"Global");
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
         if(hasEvolved)
@@ -79,119 +86,38 @@ int main(int argc, char *argv[])
         Ksl = particleCloud.momCoupleM(0).impMomSource();
         Ksl.correctBoundaryConditions();
 
+       //Force Checks
+       vector fTotal(0,0,0);
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+       reduce(fImpTotal, sumOp<vector>());
+       Info << "TotalForceExp: " << fTotal << endl;
+       Info << "TotalForceImp: " << fImpTotal << endl;
 
         #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
 
-        // get scalar source from DEM        
-        particleCloud.forceM(1).manipulateScalarField(Tsource);
-        Tsource.correctBoundaryConditions();
-
-        // solve scalar transport equation
-        fvScalarMatrix TEqn
-        (
-           fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
-         + fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
-         - fvm::laplacian(DT*voidfraction, T)
-         ==
-           Tsource
-        );
-        TEqn.relax();
-        TEqn.solve();
+        particleCloud.clockM().start(26,"Flow");
+	
+        #include "TEqn.H"
 
         if(particleCloud.solveFlow())
         {
             // Pressure-velocity PISO corrector
             {
                 // Momentum predictor
-                fvVectorMatrix UEqn
-                (
-                    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
-                  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
-//                + turbulence->divDevReff(U)
-                  + particleCloud.divVoidfractionTau(U, voidfraction)
-                  ==
-                  - fvm::Sp(Ksl/rho,U)
-                );
-
-                UEqn.relax();
-                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
-                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else if (momentumPredictor)
-                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+                 #include "UEqn.H"
 
                 // --- PISO loop
 
-                //for (int corr=0; corr<nCorr; corr++)
-                int nCorrSoph = nCorr + 5 * 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,31 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+PFLAGS+= -Dcompre
+
+EXE_INC = \
+     $(PFLAGS) \
+    -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,8 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -9,9 +12,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)
     );
 
 //========================

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"),

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();

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
@@ -83,7 +120,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 +140,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,38 @@
 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$(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

@@ -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;

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

@@ -37,10 +37,10 @@ Description
 
 #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>
 
 #define TOLERANCE_ORTHO 1e-10
 
@@ -58,11 +58,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,15 +70,13 @@ 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
@@ -93,18 +91,12 @@ inline bool spheroidGeometry(double radius, double aspectRatio,  //inputs
     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,12 +109,12 @@ 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
@@ -135,36 +127,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,22 +191,22 @@ 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
@@ -244,7 +234,7 @@ inline bool permutationDotDyadic(
 
 //--------------------------------------------------------------------
 // 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]
                                       )
@@ -260,6 +250,6 @@ inline bool doubleDotTensor333Tensor33(double tensor333[3][3][3],
 }
 
 
-}; //end of namespace
+} //end of namespace
 
 #endif

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.02";
+word compatibleLIGGGHTSversion="PFM 17.02";
+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,7 @@ Foam::cfdemCloud::cfdemCloud
     solveFlow_(true),
     verbose_(false),
     ignore_(false),
+    limitDEMForces_(false),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
     velocities_(NULL),
@@ -93,11 +98,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 +114,7 @@ Foam::cfdemCloud::cfdemCloud
     treatVoidCellsAsExplicitForce_(false),
     useDDTvoidfraction_(false),
     ddtVoidfraction_
-    (   
+    (
         IOobject
         (
             "ddtVoidfraction",
@@ -121,15 +128,7 @@ Foam::cfdemCloud::cfdemCloud
     ),
     turbulence_
     (
-        #if defined(version21) || defined(version16ext)
-            #ifdef compre
-                mesh.lookupObject<compressible::turbulenceModel>
-            #else
-                mesh.lookupObject<incompressible::turbulenceModel>
-            #endif
-        #elif defined(version15)
-            mesh.lookupObject<incompressible::RASModel>
-        #endif
+        mesh.lookupObject<turbulenceModel>
         (
             turbulenceModelType_
         )
@@ -230,12 +229,17 @@ Foam::cfdemCloud::cfdemCloud
         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()];
@@ -273,12 +277,23 @@ 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); 
+    if (!cgOK_ && cg_ > 1) FatalError<< "at least one of your models is not fit for cg !!!"<< abort(FatalError);
 }
 
 // * * * * * * * * * * * * * * * * Destructors  * * * * * * * * * * * * * * //
-Foam::cfdemCloud::~cfdemCloud()
+cfdemCloud::~cfdemCloud()
 {
     clockM().evalPar();
     clockM().normHist();
@@ -297,8 +312,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 +324,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 +344,7 @@ void Foam::cfdemCloud::giveDEMdata()
 
 // * * * * * * * * * * * * * * * protected Member Functions  * * * * * * * * * * * * * //
 
-void Foam::cfdemCloud::setNumberOfParticles(int nP)
+void cfdemCloud::setNumberOfParticles(int nP)
 {
     if(nP != numberOfParticles())
     {
@@ -337,12 +353,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 +366,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 +399,7 @@ void Foam::cfdemCloud::setParticleForceField()
     );
 }
 
-void Foam::cfdemCloud::setVectorAverages()
+void cfdemCloud::setVectorAverages()
 {
     if(verbose_) Info << "- setVectorAverage(Us,velocities_,weights_)" << endl;
     averagingM().setVectorAverage
@@ -384,12 +413,12 @@ 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)
     {
@@ -400,49 +429,61 @@ void Foam::cfdemCloud::setPos(double**& pos)
 }
 // * * * * * * * * * * * * * * * ACCESS  * * * * * * * * * * * * * //
 
-label Foam::cfdemCloud::particleCell(int index)
+label cfdemCloud::particleCell(int index)
 {
     label cellI = cellIDs()[index][0];
     return cellI;
 }
 
-vector Foam::cfdemCloud::position(int index)
+vector cfdemCloud::position(int index)
 {
     vector pos;
     for(int i=0;i<3;i++) pos[i] = positions()[index][i];
     return pos;
 }
 
-vector Foam::cfdemCloud::velocity(int index)
+vector cfdemCloud::velocity(int index)
 {
     vector vel;
     for(int i=0;i<3;i++) vel[i] = velocities()[index][i];
     return vel;
 }
 
-vector Foam::cfdemCloud::fluidVel(int index)
+vector cfdemCloud::expForce(int index)
+{
+    vector force;
+    for(int i=0;i<3;i++) force[i] = DEMForces()[index][i];
+    return force;
+}
+
+vector cfdemCloud::fluidVel(int index)
 {
     vector vel;
     for(int i=0;i<3;i++) vel[i] = fluidVels()[index][i];
     return vel;
 }
 
-const forceModel& Foam::cfdemCloud::forceM(int i)
+const forceModel& cfdemCloud::forceM(int i)
 {
     return forceModel_[i];
 }
 
-int Foam::cfdemCloud::nrForceModels()
+int cfdemCloud::nrForceModels()
 {
     return forceModels_.size();
 }
 
-scalar Foam::cfdemCloud::voidfraction(int index)
+int cfdemCloud::nrMomCoupleModels()
+{
+    return momCoupleModels_.size();
+}
+
+scalar cfdemCloud::voidfraction(int index)
 {
     return voidfractions()[index][0];
 }
 
-label Foam::cfdemCloud::liggghtsCommandModelIndex(word name)
+label cfdemCloud::liggghtsCommandModelIndex(word name)
 {
     int index=-1;
     forAll(liggghtsCommandModelList_,i)
@@ -456,12 +497,12 @@ label Foam::cfdemCloud::liggghtsCommandModelIndex(word name)
     return index;
 }
 
-std::vector<double*>* Foam::cfdemCloud::getVprobe()
+std::vector< std::vector<double*> >* cfdemCloud::getVprobe()
 {
  return probeModel_->getVprobe();
 }
 
-std::vector<double>* Foam::cfdemCloud::getSprobe()
+std::vector< std::vector<double> >* cfdemCloud::getSprobe()
 {
  return probeModel_->getSprobe();
 }
@@ -470,7 +511,7 @@ std::vector<double>* Foam::cfdemCloud::getSprobe()
 
 // * * *   write cfdemCloud internal data   * * * //
 
-bool Foam::cfdemCloud::evolve
+bool cfdemCloud::evolve
 (
     volScalarField& alpha,
     volVectorField& Us,
@@ -480,20 +521,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,21 +572,21 @@ 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;
@@ -591,6 +632,8 @@ bool Foam::cfdemCloud::evolve
             clockM().start(23,"giveDEMdata");
             giveDEMdata();
             clockM().stop("giveDEMdata");
+
+            dataExchangeM().couple(1);
         }//end dataExchangeM().couple()
 
 
@@ -607,7 +650,7 @@ bool Foam::cfdemCloud::evolve
     return doCouple;
 }
 
-bool Foam::cfdemCloud::reAllocArrays() const
+bool cfdemCloud::reAllocArrays()
 {
     if(numberOfParticlesChanged_ && !arraysReallocated_)
     {
@@ -622,7 +665,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,7 +675,7 @@ bool Foam::cfdemCloud::reAllocArrays() const
     return false;
 }
 
-bool Foam::cfdemCloud::reAllocArrays(int nP, bool forceRealloc) const
+bool cfdemCloud::reAllocArrays(int nP, bool forceRealloc)
 {
     if( (numberOfParticlesChanged_ && !arraysReallocated_) || forceRealloc)
     {
@@ -646,7 +689,7 @@ bool Foam::cfdemCloud::reAllocArrays(int nP, bool forceRealloc) const
         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(cellIDs_,-1,voidFractionM().maxCellsPerParticle(),nP);
         dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle(),nP);
         dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle(),nP);
         arraysReallocated_ = true;
@@ -687,7 +730,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 +768,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,10 @@ protected:
 
     bool ignore_;
 
+    bool limitDEMForces_;
+
+    scalar maxDEMForce_;
+
     const word modelType_;
 
     mutable double **positions_;
@@ -117,7 +118,7 @@ protected:
 
     mutable double **voidfractions_;
 
-    mutable double **cellIDs_;
+    mutable int **cellIDs_;
 
     mutable double **particleWeights_;
 
@@ -127,6 +128,8 @@ protected:
 
     int numberOfParticles_;
 
+    scalar d32_;
+
     bool numberOfParticlesChanged_;
 
     mutable bool arraysReallocated_;
@@ -137,6 +140,8 @@ protected:
 
     const wordList liggghtsCommandModelList_;
 
+    const wordList otherForceModels_;
+
     const word turbulenceModelType_;
 
     mutable scalar cg_;
@@ -155,15 +160,7 @@ protected:
 
     mutable volScalarField ddtVoidfraction_;
 
-    #if defined(version21) || defined(version16ext)
-        #ifdef compre
-            const compressible::turbulenceModel& turbulence_;
-        #else
-            const incompressible::turbulenceModel& turbulence_;
-        #endif
-    #elif defined(version15)
-        const incompressible::RASModel& turbulence_;
-    #endif
+    const turbulenceModel& turbulence_;
 
     autoPtr<forceModel>* forceModel_;
 
@@ -189,6 +186,8 @@ protected:
 
     autoPtr<liggghtsCommandModel>* liggghtsCommand_;
 
+    autoPtr<otherForceModel>* otherForceModel_;
+
 // Protected member functions
     virtual void getDEMdata();
 
@@ -229,7 +228,7 @@ public:
 
        void setPos(double **&);
 
-       word modelType(){ return modelType_; };
+       word modelType(){ return modelType_; }
 
        label particleCell(int);
 
@@ -237,12 +236,16 @@ public:
 
        vector velocity(int);
 
+       vector expForce(int);
+
        vector fluidVel(int);
 
        virtual const forceModel& forceM(int);
 
        virtual int nrForceModels();
 
+       virtual int nrMomCoupleModels();
+
        scalar voidfraction(int);
 
        label liggghtsCommandModelIndex(word);
@@ -291,9 +294,9 @@ 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;
 
@@ -305,24 +308,25 @@ public:
 
        virtual inline double d(int);
 
+       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 Foam::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 Foam::vector(0,0,0);}
 
        //Detector if SRF module is enable or not
        virtual inline bool SRFOn(){return false;}
@@ -361,24 +365,18 @@ 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 void postFlow() {}
 
-        virtual bool reAllocArrays(int nP, bool forceRealloc) const;  //force number of particles during reallocation
+        virtual bool reAllocArrays();
+
+        virtual bool reAllocArrays(int nP, bool forceRealloc);  //force number of particles during reallocation
 
 
     // IO
@@ -399,9 +397,11 @@ public:
 
         void resetArray(double**&,int,int,double resetVal=0.);
 
-        std::vector<double*>* getVprobe();
-         
-        std::vector<double>* getSprobe();
+        std::vector< std::vector<double*> >* getVprobe();
+
+        std::vector< std::vector<double> >* getSprobe();
+
+        void otherForces(volVectorField&);
 };
 
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -45,25 +45,25 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 inline void cfdemCloud::setCG(double cg) const
-{ 
+{
     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
 {
-	    return imExSplitFactor_;
-};
+    return imExSplitFactor_;
+}
 
 inline const bool& cfdemCloud::treatVoidCellsAsExplicitForce() const
 {
@@ -71,9 +71,9 @@ inline const bool& cfdemCloud::treatVoidCellsAsExplicitForce() const
 }
 
 inline const scalar& cfdemCloud::cg() const
-{ 
+{
     return cg_;
-};
+}
 
 inline const bool& cfdemCloud::ignore() const
 {
@@ -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,7 +177,7 @@ inline double ** cfdemCloud::particleWeights() const
     return particleWeights_;
 }
 
-inline label Foam::cfdemCloud::body(int index)
+inline label cfdemCloud::body(int index)
 {
     return index;
 }
@@ -197,6 +197,25 @@ inline double cfdemCloud::d(int index)
     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
 {
     return numberOfParticles_;
@@ -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)
 {
     // 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)
 {
     // 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  * * * * * * * * * * * * * //
+int cfdemCloudEnergy::nrEnergyModels()
+{
+    return energyModels_.size();
+}
+
+bool& cfdemCloudEnergy::implicitEnergyModel()
+{
+    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();
+
+       int nrEnergyModels();
+
+       inline const wordList& energyModels();
+       
+       bool& implicitEnergyModel();
+       
+       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()
+{
+    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())
     {
+        // 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,18 +216,18 @@ 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)
 {
     vector vel;
-    for(int i=0;i<3;i++) vel[i] = angularVelocities_[index][i]; 
+    for(int i=0;i<3;i++) vel[i] = angularVelocities_[index][i];
     return vel;
 }
 

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,8 +83,8 @@ 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
@@ -95,7 +95,7 @@ public:
         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);
@@ -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,7 +187,7 @@ bool cfdemCloudMS::evolve
     return false;
 }
 
-bool cfdemCloudMS::reAllocArrays() const
+bool cfdemCloudMS::reAllocArrays()
 {
     if(cfdemCloud::reAllocArrays())
     {
@@ -203,18 +196,11 @@ bool cfdemCloudMS::reAllocArrays() const
         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(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");
@@ -223,7 +209,7 @@ bool cfdemCloudMS::reAllocArrays() const
     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();
@@ -164,6 +163,8 @@ public:
 
        const forceModel& forceM(int);
 
+       bool reAllocArrays();
+
        bool evolve(volScalarField&,volVectorField&,volVectorField&);
 };
 

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

@@ -29,12 +29,14 @@ 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)
 {
     return bodies_[0][index]-1;
 }
@@ -54,31 +56,31 @@ inline double cfdemCloudMS::particleVolume(int index)
     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)
 {
     return nrigids_[0][index];
 }
 
-inline int Foam::cfdemCloudMS::nrForceModels()
+inline int cfdemCloudMS::nrForceModels()
 {
     return forceModels_.size();
 }

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

@@ -1,23 +1,2 @@
-#define version23  // currently being tested
-//#define version22  // currently being used
-//#define version21
-//#define version16ext
-//#define version15
+#define version4x
 
-// 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/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 = \

src/lagrangian/cfdemParticle/etc/bashrc

@@ -26,7 +26,8 @@
 #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_LIGGGHTS_BIN_DIR=$HOME/LIGGGHTS/LIGGGHTS-PUBLIC/src-build
+#export CFDEM_LIGGGHTS_MAKEFILE_NAME=fedora_fpic # obsolete with CMake
 #export CFDEM_LPP_DIR=$HOME/LIGGGHTS/mylpp/src
 #export CFDEM_PIZZA_DIR=$HOME/LIGGGHTS/PIZZA/gran_pizza_17Aug10/src
 #. $CFDEM_bashrc
@@ -41,7 +42,7 @@
 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 ."
+    echo "using CFDEM_LAMMPS_LIB_DIR=$CFDEM_LAMMPS_LIB_DIR defined by user."
 fi
 
 #- LIGGGHTS lib name
@@ -53,6 +54,13 @@ 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_SRC_DIR/lagrangian/cfdemParticle/etc
+else
+    echo "using CFDEM_ADD_LIBS_DIR=$CFDEM_ADD_LIBS_DIR defined by user."
+fi
+
 #-----------------------------------------------------
 # additional libraries
 

src/lagrangian/cfdemParticle/etc/cfdemSystemTest.sh

@@ -11,10 +11,15 @@ source $CFDEM_SRC_DIR/lagrangian/cfdemParticle/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"
@@ -29,8 +34,10 @@ 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"
 checkDirComment "$CFDEM_LPP_DIR" '$CFDEM_LPP_DIR' "yes"
+checkDirComment "$CFDEM_ADD_LIBS_DIR" '$CFDEM_ADD_LIBS_DIR' "yes"
 checkDirComment "$CFDEM_PIZZA_DIR" '$CFDEM_PIZZA_DIR' "no"
 checkDirComment "$CFDEM_TEST_HARNESS_PATH" '$CFDEM_TEST_HARNESS_PATH' "no"
+#checkDirComment "$C3PO_SRC_DIR" '$C3PO_SRC_DIR' "no"
 echo ""
 
 echo "library names"
@@ -61,3 +68,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

src/lagrangian/cfdemParticle/etc/cshrc

@@ -37,6 +37,13 @@
 #- export environment variables (adapt to your paths)
 #------------------------------------------------------------------------------
 
+#check if default lammps lib path should be used
+if ( ! ($?CFDEM_LAMMPS_LIB_DIR) ) then
+    setenv CFDEM_LAMMPS_LIB_DIR $CFDEM_LIGGGHTS_SRC_DIR"/../lib/"
+else
+    echo "using CFDEM_LAMMPS_LIB_DIR=$CFDEM_LAMMPS_LIB_DIR defined by user."
+endif
+
 #- LIGGGHTS lib name
 setenv CFDEM_LIGGGHTS_LIB_NAME lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME
 
@@ -46,6 +53,16 @@ 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_SRC_DIR/lagrangian/cfdemParticle/etc
+else
+    echo "using CFDEM_ADD_LIBS_DIR=$CFDEM_ADD_LIBS_DIR defined by user."
+endif
+
+#-----------------------------------------------------
+# additional libraries
+
 #- 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

src/lagrangian/cfdemParticle/etc/functions.sh

@@ -76,17 +76,17 @@ compileLib()
         i=$((${#str}-4))
         ending=${str:$i:4}
         if [[ $ending == "Comp" ]]; then
-                echo "Compiling a compressible library - so doing an rmdepall of incomp library first."
-                echo "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
 }
 
 #==================================#
@@ -285,7 +284,7 @@ cleanCFDEM()
 
             cd  $path
             echo "cleaning library $PWD"
-            rmdepall
+            wrmdep -a
             wclean    
             rm -r ./Make/linux*
             rm -r ./lnInclude
@@ -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()

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

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

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

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

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

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

src/lagrangian/cfdemParticle/subModels/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()
+{
+    delete partTemp_;
+    delete partRho_;
+
+    for (int i=0; i<speciesNames_.size();i++) delete [] concentrations_[i];
+    for (int i=0; i<speciesNames_.size();i++) delete [] changeOfSpeciesMass_[i];
+
+}
+
+// * * * * * * * * * * * * * * * 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,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
-#include "mpi.h"
+#include <mpi.h>
 #include "clockModel.H"
 #include <unistd.h>
 
@@ -50,242 +50,185 @@ defineRunTimeSelectionTable(clockModel, dictionary);
 
 void Foam::clockModel::start(int pos) const
 {
-	start(pos,"");
-	return;
+    start(pos,"");
+    return;
 }
 
-void Foam::clockModel::start(int pos,std::string ident) const
+void Foam::clockModel::start(int pos, const std::string& ident) const
 {
-    if(particleCloud_.mesh().time().value() > startTime_)
+    if (particleCloud_.mesh().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_ = 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();
     }
-	return;
+    return;
 }
 
 void Foam::clockModel::stop() const
 {
-    if(particleCloud_.mesh().time().value() > startTime_)
+    if (particleCloud_.mesh().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;
+    return;
 }
 
-void Foam::clockModel::stop(std::string ident) const
+void Foam::clockModel::stop(const std::string& ident) const
 {
-    if(particleCloud_.mesh().time().value() > startTime_)
+    if (particleCloud_.mesh().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;
+    return;
 }
 
 std::string Foam::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
 {
-	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
 {
-	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();
-	}
+        outFile.open(fileName.data(),ios_base::trunc);
+        outFile << strs.str();
+        outFile.close();
+    }
 
-	return;
+    return;
 }
 
 
@@ -294,108 +237,85 @@ void Foam::clockModel::initElems()
     //init elems
     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] = "";
+        nOfRuns_[i] = 0;
+        level_[i] = -1;
+        parent_[i] = -2;
     }
 }
 
 std::vector<int> Foam::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
 {
-	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();
+    return;
 }
 
-void Foam::clockModel::plotHist(double buffIn,std::string identifier,int numprocs,int myrank) const
+void Foam::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) globalTime_all = new double[numprocs];
+    MPI_Gather(&buffIn, 1, MPI_DOUBLE, globalTime_all, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
 
-    if(myrank==0)
-        for(int j=0;j<numprocs;j++)
+    if (myrank == 0)
+        for (int j=0; j<numprocs; j++)
             printf("%4f  ",globalTime_all[j]);
-    Info << "\t" <<identifier << endl;
 
-    particleCloud_.dataExchangeM().destroy(globalTime);
-    particleCloud_.dataExchangeM().destroy(globalTime_all);*/
+    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
@@ -406,83 +326,80 @@ void Foam::clockModel::Hist() const
     //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';
+
+    return;
 }
 
 void Foam::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("");
+    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 = static_cast<double>(SwapMem)/1024.0; //kB -> MB
+    double RssMB = static_cast<double>(RssMem)/1024.0;
 
-	//Pout << "SWAP Memory used: " << SwapMem <<"MB\n";
-	//Pout << "Rss Memory used: " << RssMem <<"MB\n";
-	
+    inFile.close();
 
-	return;
+    // 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";
+
+    return;
 }
 // * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * //
 
@@ -501,16 +418,15 @@ Foam::clockModel::clockModel
     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_)),
+    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;
 }
 
 

src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/clockModel.H

@@ -120,19 +120,19 @@ public:
 
     // Member Functions
 
-	    virtual void start(int pos) const;	                //start measurement with custom string identifier
-	    virtual void start(int pos,std::string identifier) const;   //start measurement with custom string identifier
-	    virtual void stop() const;				            //stop last started measurement
-	    virtual void stop(std::string identifier) const;    //stop last started measurement with check if identifier is equal
-	    virtual std::string eval() const;
-	    virtual void evalFile() const;
-	    virtual void evalPar() const;
-	    void initElems();
-	    std::vector<int> calcShift() const;			        //detects empty indices in vector, when times are evaluated
-	    void Hist() const;					                //calc Histogram
-	    virtual void normHist() const;				        //calc normalized Histogram
-        void plotHist(double,std::string,int,int) const;    //plot histogramm to terminal
-    	void getRAMUsage() const;
+        virtual void start(int pos) const;                  //start measurement with custom string identifier
+        virtual void start(int pos,const std::string& identifier) const;   //start measurement with custom string identifier
+        virtual void stop() const;                          //stop last started measurement
+        virtual void stop(const std::string& identifier) const;    //stop last started measurement with check if identifier is equal
+        virtual std::string eval() const;
+        virtual void evalFile() const;
+        virtual void evalPar() const;
+        void initElems();
+        std::vector<int> calcShift() const;	                //detects empty indices in vector, when times are evaluated
+        void Hist() const;                                  //calc Histogram
+        virtual void normHist() const;                      //calc normalized Histogram
+        void plotHist(double,const std::string&,int,int) const;    //plot histogramm to terminal
+        void getRAMUsage() const;
 
 };
 

src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.H

@@ -77,14 +77,14 @@ public:
 
     // Member Functions
 
-	    void start(int pos) const {};
-	    void start(int pos,std::string identifier) const {};
-	    void stop() const {};
-	    void stop(std::string identifier) const {};
-	    std::string eval() const {return "";};
-	    void evalFile() const {};
-	    void evalPar() const {};
-        void normHist() const {};
+        void start(int pos) const {}
+        void start(int pos,const std::string& identifier) const {}
+        void stop() const {}
+        void stop(const std::string& identifier) const {}
+        std::string eval() const { return ""; }
+        void evalFile() const {}
+        void evalPar() const {}
+        void normHist() const {}
 };
 
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.C

@@ -63,18 +63,17 @@ void Foam::dataExchangeModel::allocateArray
     int length
 ) const
 {
-    // this model should only be used for VTK (and File exchange model)
-    //if(particleCloud_.dataExchangeM().type() != "oneWayVTK")
-    //    FatalError<< "dataExchangeModel::allocateArray should not be used with your dataExchangeModel" << abort(FatalError);
-    
     // allocate and init double array
+    destroy(array, -1);
+    double *data = new double[width*length];
+    std::fill_n(data, width*length, initVal);
     array = new double*[length];
+
+    int n = 0;
     for (int i=0; i<length; i++)
     {
-        array[i] = new double [width];
-
-        // init array
-        for(int j=0; j<width;j++) array[i][j] = initVal;
+        array[i] = &data[n];
+        n += width;
     }
 }
 
@@ -87,17 +86,17 @@ void Foam::dataExchangeModel::allocateArray
 ) const
 {
     int len=0;
-    if(strcmp(length,"nparticles")==0) len = particleCloud_.numberOfParticles();
+    if (strcmp(length,"nparticles")==0) len = particleCloud_.numberOfParticles();
     else if (strcmp(length,"nbodies")==0) len = particleCloud_.numberOfClumps();
     else FatalError<<"call allocateArray with length, nparticles or nbodies!\n" << abort(FatalError);
     allocateArray(array,initVal,width,len);
 }
-void Foam::dataExchangeModel::destroy(double** array,int len) const
+
+void Foam::dataExchangeModel::destroy(double** array,int /*len*/) const
 {
     if (array == NULL) return;
 
-    for ( int i = 0; i < len; i++ )
-        delete [] array[i];  
+    if(array[0]) delete [] array[0];
 
     delete [] array;
 }
@@ -112,18 +111,17 @@ void Foam::dataExchangeModel::allocateArray
     int length
 ) const
 {
-    // this model should only be used for VTK (and File exchange model)
-    //if(particleCloud_.dataExchangeM().type() != "oneWayVTK")
-    //    FatalError<< "dataExchangeModel::allocateArray should not be used with your dataExchangeModel" << abort(FatalError);
-
-    // allocate and init double array
+    // allocate and init int array
+    destroy(array, -1);
+    int *data = new int[width*length];
+    std::fill_n(data, width*length, initVal);
     array = new int*[length];
+
+    int n = 0;
     for (int i=0; i<length; i++)
     {
-        array[i] = new int [width];
-
-        // init array
-        for(int j=0; j<width;j++) array[i][j] = initVal;
+        array[i] = &data[n];
+        n += width;
     }
 }
 
@@ -136,17 +134,17 @@ void Foam::dataExchangeModel::allocateArray
 ) const
 {
     int len=0;
-    if(strcmp(length,"nparticles")==0) len = particleCloud_.numberOfParticles();
+    if (strcmp(length,"nparticles")==0) len = particleCloud_.numberOfParticles();
     else if (strcmp(length,"nbodies")==0) len = particleCloud_.numberOfClumps();
     else FatalError<<"call allocateArray with length, nparticles or nbodies!\n" << abort(FatalError);
     allocateArray(array,initVal,width,len);
 }
-void Foam::dataExchangeModel::destroy(int** array,int len) const
+
+void Foam::dataExchangeModel::destroy(int** array,int /*len*/) const
 {
     if (array == NULL) return;
 
-    for ( int i = 0; i < len; i++ )
-        delete [] array[i];  
+    if(array[0]) delete [] array[0];
 
     delete [] array;
 }
@@ -161,15 +159,13 @@ void Foam::dataExchangeModel::allocateArray
     int length
 ) const
 {
-    // this model should only be used for VTK (and File exchange model)
-    //if(particleCloud_.dataExchangeM().type() != "oneWayVTK")
-    //    FatalError<< "dataExchangeModel::allocateArray should not be used with your dataExchangeModel" << abort(FatalError);
-
+    destroy(array);
     // allocate and init int array
     array = new int[length];
     for (int i=0; i<length; i++)
         array[i] = initVal;
 }
+
 void Foam::dataExchangeModel::destroy(int* array) const
 {
     if (array == NULL) return;
@@ -186,24 +182,23 @@ void Foam::dataExchangeModel::allocateArray
     int length
 ) const
 {
-    // this model should only be used for VTK (and File exchange model)
-    //if(particleCloud_.dataExchangeM().type() != "oneWayVTK")
-    //    FatalError<< "dataExchangeModel::allocateArray should not be used with your dataExchangeModel" << abort(FatalError);
-
+    destroy(array);
     // allocate and init double array
     array = new double[length];
     for (int i=0; i<length; i++)
         array[i] = initVal;
 }
+
 void Foam::dataExchangeModel::destroy(double* array) const
 {
     if (array == NULL) return;
+
     delete [] array;
 }
 //====
 
 
-bool Foam::dataExchangeModel::couple() const
+bool Foam::dataExchangeModel::couple(int i) const
 {
     bool coupleNow = false;
     if (doCoupleNow())
@@ -220,10 +215,11 @@ scalar Foam::dataExchangeModel::timeStepFraction() const
     //scalar DEMtime = DEMts_ * couplingInterval_;
     //scalar frac = ( ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() ) - (couplingStep_) * DEMtime) / DEMtime; //Chr 05.03.2013
     scalar frac = ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() - couplingStep_ * couplingTime() ) / couplingTime();
-    if (frac<1e-4) frac = 1;
+    if (frac < 1e-4) frac = 1.;
 
     return frac;
 }
+
 int Foam::dataExchangeModel::getNumberOfParticles() const
 {
     Warning << "ask for nr of particles - which is not supported for this dataExchange model" << endl;
@@ -235,6 +231,7 @@ int Foam::dataExchangeModel::getNumberOfClumps() const
     Warning << "ask for nr of clumps - which is not supported for this dataExchange model" << endl;
     return -1;
 }
+
 int Foam::dataExchangeModel::getNumberOfTypes() const
 {
     Warning << "ask for nr of types - which is not supported for this dataExchange model" << endl;

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H

@@ -69,7 +69,6 @@ protected:
         int couplingInterval_;
 
     // Protected member functions
-        void setNumberOfParticles(int) const;
 
 public:
 
@@ -116,7 +115,9 @@ public:
 
 
     // Member Function
-        inline const int& maxNumberOfParticles() const {return maxNumberOfParticles_;};
+        void setNumberOfParticles(int) const;
+
+        inline const int& maxNumberOfParticles() const { return maxNumberOfParticles_; }
 
         template <typename T>
         void getData
@@ -140,7 +141,7 @@ public:
             word type,
             int ** const& field,
             label step
-        ) const=0;
+        ) const = 0;
 
         virtual void giveData
         (
@@ -175,25 +176,25 @@ public:
         virtual void destroy(double*) const;
         //====
 
-        virtual bool couple() const;
+        virtual bool couple(int) const;
 
         virtual scalar timeStepFraction() const;
 
-        inline int couplingStep() const {return couplingStep_;};
+        inline int couplingStep() const { return couplingStep_; }
 
-        inline const scalar& DEMts() const {return DEMts_;};
+        inline const scalar& DEMts() const { return DEMts_; }
 
-        inline int couplingInterval() const {return couplingInterval_;};
+        inline int couplingInterval() const { return couplingInterval_; }
 
-        inline scalar couplingTime() const {return couplingInterval_*DEMts_;};
+        inline scalar couplingTime() const { return couplingInterval_*DEMts_; }
 
         inline scalar TSstart() const { return particleCloud_.mesh().time().startTime().value()
-                                        + (couplingStep_-1) * couplingTime();};
+                                        + (couplingStep_-1) * couplingTime(); }
 
-        inline scalar TSend() const {return particleCloud_.mesh().time().startTime().value()
-                                        + (couplingStep_) * couplingTime();};
+        inline scalar TSend() const { return particleCloud_.mesh().time().startTime().value()
+                                        + (couplingStep_) * couplingTime(); }
 
-        inline int DEMstepsTillT(scalar t) const {return (t - (particleCloud_.mesh().time().value() - couplingTime())  + SMALL) / DEMts_;};
+        inline int DEMstepsTillT(scalar t) const { return (t - (particleCloud_.mesh().time().value() - couplingTime())  + SMALL) / DEMts_; }
 
         inline void checkTSsize() const
         {
@@ -253,7 +254,7 @@ public:
 
         virtual word myType() const=0;
 
-        virtual void setCG() const { Warning << "setCG() not executed correctly!" << endl; };
+        virtual void setCG() const { Warning << "setCG() not executed correctly!" << endl; }
 };
 
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/noDataExchange/noDataExchange.H

@@ -82,7 +82,7 @@ public:
             word type,
             double ** const& field,
             label step
-        ) const {};
+        ) const {}
 
         void getData
         (
@@ -90,7 +90,7 @@ public:
             word type,
             int ** const& field,
             label step
-        ) const {};
+        ) const {}
 
         void giveData
         (
@@ -98,9 +98,9 @@ public:
             word type,
             double ** const& field,
             const char* datatype = ""
-        ) const {};
+        ) const {}
 
-        word myType() const{return typeName; };
+        word myType() const { return typeName; }
 };
 
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.C

@@ -71,10 +71,6 @@ oneWayVTK::oneWayVTK
     maxNumberOfParticles_ = readScalar(propsDict_.lookup("maxNumberOfParticles"));
     setNumberOfParticles(maxNumberOfParticles_);
 
-    // make a const char* from word
-    string HH=string(filename_);
-    charFilename_=HH.c_str();
-
     Info << "relativePath_" << relativePath_ << endl;
 }
 
@@ -94,107 +90,100 @@ void oneWayVTK::getData
     label step
 ) const
 {
-if (type == "scalar-atom")
-{
-    // get path to particle VTK files
-    char index[100];
-    sprintf(index, charFilename_, step);
-    //fileName H(particleCloud_.mesh().time().path()/".."/"DEM"/"post"/index);
-    fileName H(particleCloud_.mesh().time().path()/relativePath_/index);
-    Info << "opening file: " <<H << endl;
-
-    // set file pointer
-    string HH=string(H);
-    const char * paricleFilePath=HH.c_str();
-    ifstream* inputPtr;
-    inputPtr = new ifstream(paricleFilePath);
-
-    if (name == "radius")
+    if (type == "scalar-atom")
     {
-        // read data
-        string just_read = " ";
-        while(just_read.compare(name) != 0)  *inputPtr >> just_read;   //read until we read "name"
-        *inputPtr >> just_read;                                    // skip text for dataType
-        *inputPtr >> just_read;                                    // skip text for "1"
-        *inputPtr >> just_read;                                    // skip text for "LookUp"
-        *inputPtr >> just_read;                                    // skip text for "default"
-        for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+        // get path to particle VTK files
+        char index[100];
+        sprintf(index, filename_.c_str(), step);
+        fileName paricleFilePath(particleCloud_.mesh().time().path()/relativePath_/index);
+        Info << "opening file: " << paricleFilePath << endl;
+
+        // open file
+        ifstream input(paricleFilePath.c_str());
+
+        if (!input.is_open()) cerr << "File not found!, " << paricleFilePath << endl;
+
+        if (name == "radius")
         {
-            *inputPtr >> field[index][0];
+            // read data
+            string just_read;
+            while (just_read.compare(name) != 0)  input >> just_read; //read until we read "name"
+            input >> just_read;                                    // skip text for dataType
+            input >> just_read;                                    // skip text for "1"
+            input >> just_read;                                    // skip text for "LookUp"
+            input >> just_read;                                    // skip text for "default"
+            for (int index = 0; index<particleCloud_.numberOfParticles(); ++index)
+            {
+                input >> field[index][0];
+            }
         }
+        else
+        {
+            // read data
+            string just_read;
+            while (just_read.compare(name) != 0)  input >> just_read; //read until we read "name"
+            input >> just_read;                                    // skip text for dataType
+            for (int index = 0; index<particleCloud_.numberOfParticles(); ++index)
+            {
+                input >> field[index][0];
+            }
+        }
+
+        // close inputStream
+        input.close();
+    }
+    else if (type == "vector-atom")
+    {
+        // get path to particle VTK files
+        char index[100];
+        sprintf(index, filename_.c_str(), step);
+        fileName paricleFilePath(particleCloud_.mesh().time().path()/relativePath_/index);
+        Info << "opening file: " << paricleFilePath << endl;
+
+        // open file
+        ifstream input(paricleFilePath.c_str());
+
+        if (!input.is_open()) cerr << "File not found!, " << paricleFilePath << endl;
+
+        // read position data from VTK file
+        //NP: special case as position data has no "name" in the vtk file
+        if (name == "x")
+        {
+            int numberOfParticles;  // remove this?
+            string just_read;
+            while (just_read.compare("POINTS") != 0)  input >> just_read; // read until we read "POINTS"
+            input >> numberOfParticles;                             // this is now the number of points in the file
+            input >> just_read;                                     // skip text for dataType
+
+            // give nr of particles to cloud
+            setNumberOfParticles(numberOfParticles);
+
+            // re-allocate arrays of cloud
+            particleCloud_.reAllocArrays();
+
+            for (int index = 0; index<numberOfParticles; ++index)
+            {
+                input >> field[index][0] >> field[index][1] >> field[index][2];
+            }
+        }
+        else
+        {
+            string just_read;
+            while (just_read.compare(name) != 0)  input >> just_read;  // read until we read "name"
+            input >> just_read;                                        // skip text for dataType
+            for (int index = 0; index<particleCloud_.numberOfParticles(); ++index)
+            {
+                input >> field[index][0] >> field[index][1] >> field[index][2];
+            }
+        }
+
+        // close inputStream
+        input.close();
     }
     else
     {
-        // read data
-        string just_read = " ";
-        while(just_read.compare(name) != 0)  *inputPtr >> just_read;   //read until we read "name"
-        *inputPtr >> just_read;                                    // skip text for dataType
-        for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
-        {
-            *inputPtr >> field[index][0];
-        }
+        Info << "unknown type in getData!!!" << endl;
     }
-
-    // clean up inputStream
-    delete inputPtr;
-} else if (type == "vector-atom")
-{
-    // get path to particle VTK files
-    char index[100];
-    sprintf(index, charFilename_, step);
-    //fileName H(particleCloud_.mesh().time().path()/".."/"DEM"/"post"/index);
-    fileName H(particleCloud_.mesh().time().path()/relativePath_/index);
-    Info << "opening file: " <<H << endl;
-
-    // set file pointer
-    string HH=string(H);
-    const char * paricleFilePath=HH.c_str();
-    ifstream* inputPtr;
-    inputPtr = new ifstream(paricleFilePath);
-
-    // read position data from VTK file
-    //NP: secial case as position data has no "name" in the vtk file
-    if (name == "x")
-    {
-        int numberOfParticles;  // remove this?
-
-        string just_read = " ";
-        if(!*inputPtr) cerr << "File not found!, " << H << endl;
-        while(just_read.compare("POINTS") != 0)  *inputPtr >> just_read;   //read until we read "POINTS"
-        *inputPtr >> numberOfParticles;                           //this is now the number of points in the file
-        *inputPtr >> just_read;                                    // skip text for dataType
-
-        // give nr of particles to cloud
-        setNumberOfParticles(numberOfParticles);
-
-        // re-allocate arrays of cloud
-        particleCloud_.reAllocArrays();
-
-        for(int index = 0;index <  numberOfParticles; ++index)
-        {
-            *inputPtr >> field[index][0] >> field[index][1] >> field[index][2];
-        }
-    }
-    else
-    {
-        string just_read = " ";
-        while(just_read.compare(name) != 0)  *inputPtr >> just_read;   //read until we read "name"
-        *inputPtr >> just_read;                                    // skip text for dataType
-        for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
-        {
-            *inputPtr >> field[index][0] >> field[index][1] >> field[index][2];
-        }
-    }
-
-    // clean up inputStream
-    delete inputPtr;
-}
-else
-{
-    Info << "unknown type in getData!!!" << endl;
-}
-
-
 }
 
 void oneWayVTK::giveData

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/oneWayVTK/oneWayVTK.H

@@ -64,8 +64,6 @@ private:
 
     fileName relativePath_;
 
-    const char* charFilename_;
-
 public:
 
     //- Runtime type information
@@ -102,7 +100,7 @@ public:
             int ** const& field,
             label step
         ) const
-        {};
+        {}
 
         void giveData
         (
@@ -112,7 +110,7 @@ public:
             const char* datatype = ""
         ) const;
 
-        word myType() const{return typeName; };
+        word myType() const{ return typeName; }
 };
 
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.C

@@ -80,25 +80,12 @@ twoWayFiles::~twoWayFiles()
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-const char* twoWayFiles::wordToChar(word& inWord) const
-{
-    string HH = string(inWord);
-    return HH.c_str();
-}
-
-const char* twoWayFiles::fileNameToChar(fileName& inWord) const
-{
-    string HH = string(inWord);
-    return HH.c_str();
-}
-
 fileName twoWayFiles::getFilePath(word& name, bool in) const
 {
-    const char* charName = wordToChar(name);
     char timeStep[40];
 
     // file touched by DEM
-    strcpy(timeStep, charName);
+    strcpy(timeStep, name.c_str());
     strcat(timeStep,"1");
     fileName particleFilePathOld(particleCloud_.mesh().time().path()/"couplingFiles"/timeStep);
 
@@ -107,7 +94,7 @@ fileName twoWayFiles::getFilePath(word& name, bool in) const
     {
         Info << "wait for file " << particleFilePathOld << endl;
         struct stat st;
-        while (stat(fileNameToChar(particleFilePathOld),&st)) sleep(0.03);
+        while (stat(particleFilePathOld.c_str(),&st)) sleep(0.03);
     }
 
     return particleFilePathOld;
@@ -115,16 +102,15 @@ fileName twoWayFiles::getFilePath(word& name, bool in) const
 
 void twoWayFiles::renameFilePath(fileName& particleFilePathOld,word& name) const
 {
-    const char* charName = wordToChar(name);
     char timeStep[40];
 
     // file touched by CFD
-    strcpy(timeStep, charName);
+    strcpy(timeStep, name.c_str());
     strcat(timeStep,"0");
     fileName particleFilePath(particleCloud_.mesh().time().path()/"couplingFiles"/timeStep);
 
     // rename old file
-    rename(fileNameToChar(particleFilePathOld),fileNameToChar(particleFilePath));
+    rename(particleFilePathOld.c_str(),particleFilePath.c_str());
 }
 
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayFiles/twoWayFiles.H

@@ -43,7 +43,7 @@ SourceFiles
 
 #include "dataExchangeModel.H"
 #include "OFstream.H"
-#include "sys/stat.h"
+#include <sys/stat.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,10 +65,6 @@ private:
 
   // private member functions
 
-    const char* wordToChar(word&) const;
-
-    const char* fileNameToChar(fileName&) const;
-
     fileName getFilePath(word&, bool) const;
 
     void renameFilePath(fileName&,word&) const;
@@ -108,7 +104,7 @@ public:
             word type,
             int ** const& field,
             label step
-        ) const {FatalError<<"function not implemented !!! twoWayFiles::getData!!!\n" << abort(FatalError);};
+        ) const { FatalError<<"function not implemented !!! twoWayFiles::getData!!!\n" << abort(FatalError); }
 
         void giveData
         (
@@ -118,7 +114,7 @@ public:
             const char* datatype = ""
         ) const;
 
-        word myType() const{return typeName; };
+        word myType() const{ return typeName; }
 };
 
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C

@@ -63,45 +63,20 @@ twoWayMPI::twoWayMPI
 )
 :
     dataExchangeModel(dict,sm),
-    propsDict_(dict.subDict(typeName + "Props"))
+    propsDict_(dict.subDict(typeName + "Props")),
+    lmp(NULL)
 {
-    // set max nr of particles from dict
-    //Info << "twoWayMPI.C- this should no longer be needed" << endl;
-    //maxNumberOfParticles_ = readScalar(propsDict_.lookup("maxNumberOfParticles"));
-
-
     Info<<"Starting up LIGGGHTS for first time execution"<<endl;
 
-    MPI_Comm_rank(MPI_COMM_WORLD,&me);
-    MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
+    MPI_Comm_dup(MPI_COMM_WORLD, &comm_liggghts);
 
-    if (me < nprocs) liggghts = 1;
-    else liggghts = MPI_UNDEFINED;
-
-    MPI_Comm_split(MPI_COMM_WORLD,liggghts,0,&comm_liggghts);
+    // read path from dictionary
+    const fileName liggghtsPath(propsDict_.lookup("liggghtsPath"));
 
     // open LIGGGHTS input script
-    char *liggghtsPathChar = new char[256];
-    int n = 0;
-    if (me == 0)
-    {
-      // read path from dictionary
-      const fileName liggghtsPath(propsDict_.lookup("liggghtsPath"));
-      strcpy(liggghtsPathChar, liggghtsPath.c_str());
-      n = strlen(liggghtsPathChar) + 1;
-
-      Info<<"Executing input script '"<< liggghtsPath.c_str() <<"'"<<endl;
-    }
-
-    if (liggghts == 1) lmp = new LAMMPS_NS::LAMMPS(0,NULL,comm_liggghts);
-
-    MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
-    if (n > 0) {
-        MPI_Bcast(liggghtsPathChar,n,MPI_CHAR,0,MPI_COMM_WORLD);
-        if (liggghts == 1) lmp->input->file(liggghtsPathChar);
-    }
-
-    delete [] liggghtsPathChar;
+    Info<<"Executing input script '"<< liggghtsPath.c_str() <<"'"<<endl;
+    lmp = new LAMMPS_NS::LAMMPS(0,NULL,comm_liggghts);
+    lmp->input->file(liggghtsPath.c_str());
 
     // get DEM time step size
     DEMts_ = lmp->update->dt;
@@ -112,28 +87,20 @@ twoWayMPI::twoWayMPI
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 twoWayMPI::~twoWayMPI()
-{}
-
-// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-char* twoWayMPI::wordToChar(word& inWord) const
 {
-    string HH = string(inWord);
-    return const_cast<char*>(HH.c_str());
+    delete lmp;
 }
 
-
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 void twoWayMPI::getData
 (
     word name,
     word type,
     double ** const& field,
-    label step
+    label /*step*/
 ) const
 {
-    char* charName = wordToChar(name);
-    char* charType = wordToChar(type);
-    data_liggghts_to_of(charName,charType, lmp, (void*&) field,"double");
+    data_liggghts_to_of(name.c_str(), type.c_str(), lmp, (void*&) field, "double");
 }
 
 void twoWayMPI::getData
@@ -141,12 +108,10 @@ void twoWayMPI::getData
     word name,
     word type,
     int ** const& field,
-    label step
+    label /*step*/
 ) const
 {
-    char* charName = wordToChar(name);
-    char* charType = wordToChar(type);
-    data_liggghts_to_of(charName,charType, lmp, (void*&) field,"int");
+    data_liggghts_to_of(name.c_str(), type.c_str(), lmp, (void*&) field, "int");
 }
 
 void twoWayMPI::giveData
@@ -157,11 +122,9 @@ void twoWayMPI::giveData
     const char* datatype
 ) const
 {
-    char* charName = wordToChar(name);
-    char* charType = wordToChar(type);
-    char* charDatatype= const_cast<char*> (datatype);
-    data_of_to_liggghts(charName,charType, lmp, (void*) field,charDatatype);
+    data_of_to_liggghts(name.c_str(), type.c_str(), lmp, (void*)field, datatype);
 }
+
 //============
 // double **
 void Foam::twoWayMPI::allocateArray
@@ -172,8 +135,7 @@ void Foam::twoWayMPI::allocateArray
     int length
 ) const
 {
-    //if(length==-1) then LIGGGHTS uses own length data
-    allocate_external_double(array, width,length,initVal,lmp);
+    allocate_external_double(array, width, length, initVal, lmp);
 }
 
 void Foam::twoWayMPI::allocateArray
@@ -184,20 +146,18 @@ void Foam::twoWayMPI::allocateArray
     const char* length
 ) const
 {
-    //if(length==-1) then LIGGGHTS uses own length data
-    char* charLength= const_cast<char*> (length);
-    allocate_external_double(array, width,charLength,initVal,lmp);
+    allocate_external_double(array, width, length, initVal, lmp);
 }
-void Foam::twoWayMPI::destroy(double** array,int len) const
+
+void Foam::twoWayMPI::destroy(double** array,int /*len*/) const
 {
     if (array == NULL) return;
 
-    //for ( int i = 0; i < len; i++ ) // does not work
-    for ( int i = 0; i < 1; i++ )
-        free(array[i]); 
+    if (array[0]) free(array[0]);
 
     free(array);
 }
+
 //============
 // int **
 void Foam::twoWayMPI::allocateArray
@@ -208,8 +168,7 @@ void Foam::twoWayMPI::allocateArray
     int length
 ) const
 {
-    //if(length==-1) then LIGGGHTS uses own length data
-    allocate_external_int(array, width,length,initVal,lmp);
+    allocate_external_int(array, width, length, initVal, lmp);
 }
 
 void Foam::twoWayMPI::allocateArray
@@ -220,17 +179,14 @@ void Foam::twoWayMPI::allocateArray
     const char* length
 ) const
 {
-    //if(length==-1) then LIGGGHTS uses own length data
-    char* charLength= const_cast<char*> (length);
-    allocate_external_int(array, width,charLength,initVal,lmp);
+    allocate_external_int(array, width, length, initVal, lmp);
 }
-void Foam::twoWayMPI::destroy(int** array,int len) const
+
+void Foam::twoWayMPI::destroy(int** array,int /*len*/) const
 {
     if (array == NULL) return;
 
-    //for ( int i = 0; i < len; i++ ) // does not work
-    for ( int i = 0; i < 1; i++ )
-        free(array[i]); 
+    if (array[0]) free(array[0]);
 
     free(array);
 }
@@ -250,17 +206,15 @@ void Foam::twoWayMPI::destroy(double* array) const
 }
 //============
 
-bool Foam::twoWayMPI::couple() const
+bool Foam::twoWayMPI::couple(int i) const
 {
     bool coupleNow = false;
-    if (doCoupleNow())
+    if (i==0)
     {
         couplingStep_++;
         coupleNow = true;
 
         // start liggghts
-        if (liggghts == 1)
-        {
             // run commands from liggghtsCommands dict
             Info<<"Starting up LIGGGHTS" << endl;
             particleCloud_.clockM().start(3,"LIGGGHTS");
@@ -271,7 +225,6 @@ bool Foam::twoWayMPI::couple() const
             DynamicList<scalar> interruptTimes(0);
             DynamicList<int> DEMstepsToInterrupt(0);
             DynamicList<int> lcModel(0);
-            scalar interruptTime = -1;
 
             forAll(particleCloud_.liggghtsCommandModelList(),i)
             {
@@ -303,11 +256,11 @@ bool Foam::twoWayMPI::couple() const
                         ind = len-i-1;
                         if(ind>0)
                             DEMstepsToInterrupt[ind] -= DEMstepsToInterrupt[ind-1];
-                    }                    
+                    }
 
-                    Info << "Foam::twoWayMPI::couple(): interruptTimes=" << interruptTimes << endl;
-                    Info << "Foam::twoWayMPI::couple(): DEMstepsToInterrupt=" << DEMstepsToInterrupt << endl;
-                    Info << "Foam::twoWayMPI::couple(): lcModel=" << lcModel << endl;
+                    Info << "Foam::twoWayMPI::couple(i): interruptTimes=" << interruptTimes << endl;
+                    Info << "Foam::twoWayMPI::couple(i): DEMstepsToInterrupt=" << DEMstepsToInterrupt << endl;
+                    Info << "Foam::twoWayMPI::couple(i): lcModel=" << lcModel << endl;
                 }
 
                 if(particleCloud_.liggghtsCommand()[i]().type()=="runLiggghts")
@@ -322,16 +275,16 @@ bool Foam::twoWayMPI::couple() const
                 //    sort interrupt list within this run period
                 //    keep track of corresponding liggghtsCommand
                 int DEMstepsRun(0);
-                
+
                 forAll(interruptTimes,j)
-                {                  
+                {
                     // set run command till interrupt
-                    DEMstepsRun += DEMstepsToInterrupt[j];          
+                    DEMstepsRun += DEMstepsToInterrupt[j];
                     particleCloud_.liggghtsCommand()[runComNr]().set(DEMstepsToInterrupt[j]);
                     const char* command = particleCloud_.liggghtsCommand()[runComNr]().command(0);
                     Info << "Executing run command: '"<< command <<"'"<< endl;
                     lmp->input->one(command);
-    
+
                     // run liggghts command with exact timing
                     command = particleCloud_.liggghtsCommand()[lcModel[j]]().command(0);
                     Info << "Executing command: '"<< command <<"'"<< endl;
@@ -386,7 +339,6 @@ bool Foam::twoWayMPI::couple() const
 
             particleCloud_.clockM().stop("LIGGGHTS");
             Info<<"LIGGGHTS finished"<<endl;
-        }
 
         // give nr of particles to cloud
         double newNpart = liggghts_get_maxtag(lmp);
@@ -414,7 +366,7 @@ int Foam::twoWayMPI::getNumberOfClumps() const
     #endif
 
     Warning << "liggghts_get_maxtag_ms(lmp) is not available here!" << endl;
-    return -1;       
+    return -1;
 }
 
 int Foam::twoWayMPI::getNumberOfTypes() const
@@ -433,7 +385,7 @@ double* Foam::twoWayMPI::getTypeVol() const
     #endif
 
     Warning << "liggghts_get_vclump_ms(lmp) is not available here!" << endl;
-    return NULL;       
+    return NULL;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H

@@ -44,7 +44,7 @@ SourceFiles
 #include "dataExchangeModel.H"
 #include "liggghtsCommandModel.H"
 #include "OFstream.H"
-#include "sys/stat.h"
+#include <sys/stat.h>
 #include "pair.h"
 #include "force.h"
 #include "forceModel.H"
@@ -54,7 +54,7 @@ SourceFiles
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
-#include "mpi.h"
+#include <mpi.h>
 #include <lammps.h>         // these are LAMMPS include files
 #include <input.h>
 #include <atom.h>
@@ -81,21 +81,12 @@ private:
 
   // private data
     dictionary propsDict_;
-
-    // variables
-    int me;
-
-    int nprocs;
-
-    int liggghts;
-
     MPI_Comm comm_liggghts;
 
-    LAMMPS_NS::LAMMPS *lmp;
-
   // private member functions
 
-    char* wordToChar(word&) const;
+protected:
+    LAMMPS_NS::LAMMPS *lmp;
 
 public:
 
@@ -145,12 +136,12 @@ public:
         //============
         // double **
         void allocateArray(double**&, double, int, int) const;
-        void allocateArray(double**&, double, int,const char* ="nparticles") const;
+        void allocateArray(double**&, double, int, const char* = "nparticles") const;
         void destroy(double**,int) const;
         //============
         // int **
         void allocateArray(int**&, int, int, int) const;
-        void allocateArray(int**&, int, int,const char* ="nparticles") const;
+        void allocateArray(int**&, int, int, const char* = "nparticles") const;
         void destroy(int**,int) const;
         //==============
         // double *
@@ -160,16 +151,16 @@ public:
         void destroy(int*) const;
         //==============
 
-        bool couple() const;
+        bool couple(int) const;
 
         int getNumberOfParticles() const;
         int getNumberOfClumps() const;
         int getNumberOfTypes() const;
         double* getTypeVol() const;
 
-        word myType() const{return typeName; };
+        word myType() const { return typeName; }
 
-        void setCG() const { particleCloud_.setCG(lmp->force->cg()); };
+        void setCG() const { particleCloud_.setCG(lmp->force->cg()); }
 };
 
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic

@@ -6,7 +6,7 @@ SHELL = /bin/sh
 
 CC =		mpic++
 CCFLAGS =	-O2 -fPIC \
-		 -funroll-loops -fstrict-aliasing -Wall -W -Wno-uninitialized
+		 -funroll-loops -fstrict-aliasing -Wall -Wextra -Wno-unused-result -Wno-literal-suffix
 DEPFLAGS =	-M
 LINK =		mpic++
 LINKFLAGS =	-O -fPIC

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic_debug

@@ -6,7 +6,7 @@ SHELL = /bin/sh
 
 CC =		mpic++
 CCFLAGS =	-g -O2 -fPIC \
-		 -funroll-loops -fstrict-aliasing -Wall -W -Wno-uninitialized
+		 -funroll-loops -fstrict-aliasing -Wall -Wextra -Wno-unused-result -Wno-literal-suffix
 DEPFLAGS =	-M
 LINK =		mpic++
 LINKFLAGS =	-O -fPIC

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/error.cpp

@@ -1,6 +1,6 @@
-#include "mpi.h"
+#include <mpi.h>
 #include "stdlib.h"
-#include "stdio.h"
+#include <stdio.h>
 #include "error.h"
 
 /* ---------------------------------------------------------------------- */

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/error.h

@@ -1,7 +1,7 @@
 #ifndef ERROR_H
 #define ERROR_H
 
-#include "mpi.h"
+#include <mpi.h>
 
 class Error {
  public:

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/files.cpp

@@ -1,4 +1,4 @@
-#include "stdio.h"
+#include <stdio.h>
 #include "string.h"
 #include "files.h"
 
@@ -45,4 +45,5 @@ char **extract(char *file, char *header, int n, char **lines)
   }
 
   fclose(fp);
+  return lines;
 }

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/irregular.cpp

@@ -1,4 +1,4 @@
-#include "stdio.h"
+#include <stdio.h>
 #include "stdlib.h"
 #include "string.h"
 #include "irregular.h"

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/irregular.h

@@ -1,7 +1,7 @@
 #ifndef IRREGULAR_H
 #define IRREGULAR_H
 
-#include "mpi.h"
+#include <mpi.h>
 
 class Irregular {
  public:

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/lammps_data_write.h

@@ -2,7 +2,7 @@
 #define LAMMPS_DATA_WRITE_H
 
 #include "send2one.h"
-#include "stdio.h"
+#include <stdio.h>
 
 class LAMMPSDataWrite : public Send2One {
  public:

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/many2many.cpp

@@ -1,5 +1,5 @@
-#include "mpi.h"
-#include "stdio.h"
+#include <mpi.h>
+#include <stdio.h>
 #include "many2many.h"
 #include "irregular.h"
 #include "memory.h"
@@ -268,7 +268,6 @@ void Many2Many::exchange(int *src, int *dest)
   // user src_off and dest_off to pack/unpack data
 
   if (irregular) {
-    int nrecv = irregular->size(sizeof(int)) / sizeof(int);
     for (i = 0; i < nsrc_off; i++) src_iwork[i] = src[src_off[i]];
     irregular->exchange((char *) src_iwork, (char *) dest_iwork);
     for (i = 0; i < ndest_off; i++) dest[dest_off[i]] = dest_iwork[i];
@@ -293,7 +292,6 @@ void Many2Many::exchange(double *src, double *dest)
   // user src_off and dest_off to pack/unpack data
 
   if (irregular) {
-    int nrecv = irregular->size(sizeof(double)) / sizeof(double);
     for (i = 0; i < nsrc_off; i++) src_dwork[i] = src[src_off[i]];
     irregular->exchange((char *) src_dwork, (char *) dest_dwork);
     for (i = 0; i < ndest_off; i++) dest[dest_off[i]] = dest_dwork[i];

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/many2many.h

@@ -1,7 +1,7 @@
 #ifndef MANY2MANY_H
 #define MANY2MANY_H
 
-#include "mpi.h"
+#include <mpi.h>
 
 class Many2Many {
  public:

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/many2one.cpp

@@ -1,5 +1,5 @@
-#include "mpi.h"
-#include "stdio.h"
+#include <mpi.h>
+#include <stdio.h>
 #include "stdlib.h"
 #include "many2one.h"
 #include "memory.h"

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/many2one.h

@@ -1,7 +1,7 @@
 #ifndef MANY2ONE_H
 #define MANY2ONE_H
 
-#include "mpi.h"
+#include <mpi.h>
 
 class Many2One {
  public:

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/memory.cpp

@@ -1,6 +1,6 @@
-#include "mpi.h"
+#include <mpi.h>
 #include "stdlib.h"
-#include "stdio.h"
+#include <stdio.h>
 #include "memory.h"
 #include "error.h"
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/memory.h

@@ -1,7 +1,7 @@
 #ifndef MEMORY_H
 #define MEMORY_H
 
-#include "mpi.h"
+#include <mpi.h>
 
 class Memory {
  public: