release on 2014-05-13_08-58-58

README

@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM. Note: this code is not part of OpenFOAM (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+
+CFDEM coupling provides an open source parallel coupled CFD-DEM framework 
+combining the strengths of LIGGGHTS DEM code and the Open Source 
+CFD package OpenFOAM(R)(*). The CFDEMcoupling toolbox allows to expand 
+standard CFD solvers of OpenFOAM(R)(*) to include a coupling to the DEM 
+code LIGGGHTS. In this toolbox the particle representation within the 
+CFD solver is organized by "cloud" classes. Key functionalities are organised 
+in sub-models (e.g. force models, data exchange models, etc.) which can easily 
+be selected and combined by dictionary settings.
+
+The coupled solvers run fully parallel on distributed-memory clusters. 
+
+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
+
+The file structure:
+
+- "src" directory including the source files of the coupling toolbox and models
+- "applications" directory including the solver files for coupled CFD-DEM simulations
+- "doc" directory including the documentation of CFDEMcoupling
+- "tutorials" directory including basic tutorial cases showing the functionality
+
+
+
+Details on installation are given on the "www.cfdem.com"
+
+The functionality of this CFD-DEM framwork is described via "tutorial cases" showing 
+how to use different solvers and models.
+
+CFDEMcoupling stands for Computational Fluid Dynamics (CFD) -
+Discrete Element Method (DEM) coupling.
+
+CFDEMcoupling is an open-source code, distributed freely under the terms of the 
+GNU Public License (GPL).
+
+Core development of CFDEMcoupling is done by 
+Christoph Goniva and Christoph Kloss, both at DCS Computing GmbH, 2012
+
+
+\*---------------------------------------------------------------------------*/
+(*) "OpenFOAM(R)"_of is a registered trade mark of the ESI Group.
+This offering is not affiliated, approved or endorsed by ESI Group, 
+the producer of the OpenFOAM® software and owner of the OpenFOAM® trade mark.
+\*---------------------------------------------------------------------------*/

applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C

@@ -83,95 +83,103 @@ int main(int argc, char *argv[])
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
-        // Pressure-velocity PISO corrector
+
+        if(particleCloud.solveFlow())
         {
-            // Momentum predictor
+            // Pressure-velocity PISO corrector
-            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)
-            );
-
-            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++)
             {
-                volScalarField rUA = 1.0/UEqn.A();
+                // Momentum predictor
-
+                fvVectorMatrix UEqn
-                surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+                (
-                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+                    fvm::ddt(voidfraction,U) + fvm::Sp(fvc::ddt(voidfraction),U)
-
+                  + fvm::div(phi,U) + fvm::Sp(fvc::div(phi),U)
-                U = rUA*UEqn.H();
+//                + turbulence->divDevReff(U)
-
+                  + particleCloud.divVoidfractionTau(U, voidfraction)
-                phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() );
+                  ==
-                     //+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
+                  - fvm::Sp(Ksl/rho,U)
-                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();
-                    }
-
-                } // end non-orthogonal corrector loop
-
-                #include "continuityErrorPhiPU.H"
 
                 if (modelType=="B")
-                    U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
+                 UEqn == - fvc::grad(p) + Ksl/rho*Us;
                 else
-                    U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
+                    UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
 
-                U.correctBoundaryConditions();
+                UEqn.relax();
 
-            } // end piso loop
+                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++)
+                {
+                    volScalarField rUA = 1.0/UEqn.A();
+
+                    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+                    volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+
+                    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++)
+                    {
+                        // 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();
+                        }
+
+                    } // 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
+            }
+
+            turbulence->correct();
+        }// end solveFlow
+        else
+        {
+            Info << "skipping flow solution." << endl;
         }
 
-        turbulence->correct();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

doc/forceModel_GidaspowDrag.html

@@ -24,6 +24,7 @@ GidaspowDragProps
     voidfractionFieldName "voidfraction";
     phi "scalar";
     interpolation;
+    implDEM;
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
@@ -36,6 +37,7 @@ GidaspowDragProps
 
 <LI><I>interpolation</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
 
+<I>implDEM</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 
 </UL>
 <P><B>Examples:</B>

doc/forceModel_GidaspowDrag.txt

@@ -22,6 +22,7 @@ GidaspowDragProps
     voidfractionFieldName "voidfraction";
     phi "scalar";
     interpolation;
+    implDEM;
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
@@ -29,6 +30,7 @@ GidaspowDragProps
 {voidfraction} = name of the finite volume voidfraction field :l
 {phi} = drag correction factor (in doubt 1) :l
 {interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
+{implDEM} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 :ule
 
 [Examples:]

doc/forceModel_KochHillDrag.html

@@ -23,6 +23,7 @@ KochHillDragProps
     densityFieldName "density";
     voidfractionFieldName "voidfraction"; 
     interpolation;
+    implDEM;
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
@@ -33,8 +34,7 @@ KochHillDragProps
 
 <LI><I>interpolation</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
 
-<LI><I>implDEM</I> = flag to use implicit formulation of drag on DEM side (normally off) 
+<I>implDEM</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
-
 
 </UL>
 <P><B>Examples:</B>

doc/forceModel_KochHillDrag.txt

@@ -21,13 +21,14 @@ KochHillDragProps
     densityFieldName "density";
     voidfractionFieldName "voidfraction"; 
     interpolation;
+    implDEM;
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
 {density} = name of the finite volume gravity field :l
 {voidfraction} = name of the finite volume voidfraction field :l
 {interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
-{implDEM} = flag to use implicit formulation of drag on DEM side (normally off) :l
+{implDEM} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 :ule
 
 [Examples:]

doc/momCoupleModel.html

@@ -29,6 +29,10 @@
 </PRE>
 <P>Note: This examples list might not be complete - please look for other models (momCoupleModel_XY) in this documentation.
 </P>
+<P>Forces can be coupled in an implicit way to the fluid solver (i.e., when solving the Navier-Stokes equations, the fluid velocity at the new time will be considered for the coupling force). This implicit coupling is typically done for the drag forces (look for "impForces()" in the implementation of the drag model). Implicit coupling is more stable (especially important for dense flows), but conflicts Newton's third law. Explicit forces are imposed on the flow solver in an explicit fashion (look for "expForces()" in the implementation of the drag model), which is less stable, but does not conflict Newton's third law.
+</P>
+<P>Note that the variable "imExSplitFactor" can be set in the couplingProperties in order to treat implicitly defined forces (in the implementation of the force model) as explicit ones. "imExSplitFactor 1.0;" is set by default, meaning that all implicit forces will be considered implicitly, whereas "imExSplitFactor 0.0;" would mean that implicitly defined forces will be treated in an explicit fashion.
+</P>
 <P><B>Description:</B>
 </P>
 <P>The momCoupleModel is the base class for momentum exchange between DEM and CFD simulation.

doc/momCoupleModel.txt

@@ -27,6 +27,10 @@ momCoupleModels
 
 Note: This examples list might not be complete - please look for other models (momCoupleModel_XY) in this documentation.
 
+Forces can be coupled in an implicit way to the fluid solver (i.e., when solving the Navier-Stokes equations, the fluid velocity at the new time will be considered for the coupling force). This implicit coupling is typically done for the drag forces (look for "impForces()" in the implementation of the drag model). Implicit coupling is more stable (especially important for dense flows), but conflicts Newton's third law. Explicit forces are imposed on the flow solver in an explicit fashion (look for "expForces()" in the implementation of the drag model), which is less stable, but does not conflict Newton's third law.
+
+Note that the variable "imExSplitFactor" can be set in the couplingProperties in order to treat implicitly defined forces (in the implementation of the force model) as explicit ones. "imExSplitFactor 1.0;" is set by default, meaning that all implicit forces will be considered implicitly, whereas "imExSplitFactor 0.0;" would mean that implicitly defined forces will be treated in an explicit fashion.
+
 [Description:]
 
 The momCoupleModel is the base class for momentum exchange between DEM and CFD simulation.

src/lagrangian/cfdemParticle/Make/files

@@ -6,7 +6,6 @@ voidFractionModels = subModels/voidFractionModel
 locateModels = subModels/locateModel
 meshMotionModels = subModels/meshMotionModel
 momCoupleModels = subModels/momCoupleModel
-regionModels = subModels/regionModel
 dataExchangeModels = subModels/dataExchangeModel
 averagingModels = subModels/averagingModel
 clockModels = subModels/clockModel
@@ -15,6 +14,7 @@ smoothingModels = subModels/smoothingModel
 probeModels = subModels/probeModel
 
 $(cfdemCloud)/cfdemCloud.C
+derived/cfdemCloudBiDisperse/cfdemCloudBiDisperse.C
 derived/cfdemCloudIB/cfdemCloudIB.C
 derived/cfdemCloudMS/cfdemCloudMS.C
 
@@ -23,36 +23,59 @@ $(forceModels)/forceModel/newForceModel.C
 $(forceModels)/noDrag/noDrag.C
 $(forceModels)/checkCouplingInterval/checkCouplingInterval.C
 $(forceModels)/DiFeliceDrag/DiFeliceDrag.C
+$(forceModels)/DiFeliceDragNLift/DiFeliceDragNLift.C
 $(forceModels)/GidaspowDrag/GidaspowDrag.C
 $(forceModels)/SchillerNaumannDrag/SchillerNaumannDrag.C
 $(forceModels)/Archimedes/Archimedes.C
 $(forceModels)/ArchimedesIB/ArchimedesIB.C
 $(forceModels)/interface/interface.C
 $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C
+$(forceModels)/interfaceParticleProbe/interfaceParticleProbe.C
+$(forceModels)/fieldStore/fieldStore.C
+$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
+$(forceModels)/fieldBound/fieldBound.C
+$(forceModels)/volWeightedAverage/volWeightedAverage.C
+$(forceModels)/totalMomentumExchange/totalMomentumExchange.C
 $(forceModels)/KochHillDrag/KochHillDrag.C
 $(forceModels)/KochHillRWDrag/KochHillRWDrag.C
+$(forceModels)/BeetstraDrag/multiphaseFlowBasic/multiphaseFlowBasic.C
+$(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/LaEuScalarLiquid/LaEuScalarLiquid.C
 $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
+$(forceModels)/LaEuScalarDust/LaEuScalarDust.C
 $(forceModels)/virtualMassForce/virtualMassForce.C
 $(forceModels)/gradPForce/gradPForce.C
+$(forceModels)/gradULiftForce/gradULiftForce.C
+$(forceModels)/HollowayDrag/HollowayDrag.C
 $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C
+$(forceModels)/melting/melting.C
+$(forceModels)/KochHillDragNLift/KochHillDragNLift.C
+$(forceModels)/stokesSpheroidDrag/stokesSpheroidDrag.C
+$(forceModels)/solidsPressureForce/solidsPressureForce.C
+$(forceModels)/periodicPressure/periodicPressure.C
+$(forceModels)/periodicPressureControl/periodicPressureControl.C
+$(forceModels)/averageSlipVel/averageSlipVel.C
 $(forceModels)/particleCellVolume/particleCellVolume.C
-$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
-$(forceModels)/volWeightedAverage/volWeightedAverage.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
 $(forceModelsMS)/DiFeliceDragMS/DiFeliceDragMS.C
+$(forceModelsMS)/GidaspowDragMS/GidaspowDragMS.C
+$(forceModelsMS)/noDragMS/noDragMS.C
 
 $(probeModels)/probeModel/probeModel.C
 $(probeModels)/probeModel/newProbeModel.C
 $(probeModels)/noProbe/noProbe.C
 $(probeModels)/particleProbe/particleProbe.C
+/*$(probeModels)/interfaceParticleProbe/interfaceParticleProbe.C*/
 
 $(IOModels)/IOModel/IOModel.C
 $(IOModels)/IOModel/newIOModel.C
 $(IOModels)/noIO/noIO.C
 $(IOModels)/basicIO/basicIO.C
+$(IOModels)/tempIO/tempIO.C
+$(IOModels)/colorIO/colorIO.C
 $(IOModels)/trackIO/trackIO.C
 $(IOModels)/sophIO/sophIO.C
 
@@ -60,34 +83,36 @@ $(voidFractionModels)/voidFractionModel/voidFractionModel.C
 $(voidFractionModels)/voidFractionModel/newVoidFractionModel.C
 $(voidFractionModels)/centreVoidFraction/centreVoidFraction.C
 $(voidFractionModels)/dividedVoidFraction/dividedVoidFraction.C
+$(voidFractionModels)/dividedVoidFractionBiDi/dividedVoidFractionBiDi.C
 $(voidFractionModels)/dividedVoidFractionMS/dividedVoidFractionMS.C
 $(voidFractionModels)/bigParticleVoidFraction/bigParticleVoidFraction.C
 $(voidFractionModels)/GaussVoidFraction/GaussVoidFraction.C
 $(voidFractionModels)/IBVoidFraction/IBVoidFraction.C
+$(voidFractionModels)/weightedNeigbhorsVoidFraction/weightedNeigbhorsVoidFraction.C
 
 $(locateModels)/locateModel/locateModel.C
 $(locateModels)/locateModel/newLocateModel.C
 $(locateModels)/standardSearch/standardSearch.C
 $(locateModels)/engineSearch/engineSearch.C
+$(locateModels)/engineSearchMany2Many/engineSearchMany2Many.C
 $(locateModels)/turboEngineSearch/turboEngineSearch.C
 $(locateModels)/turboEngineSearchM2M/turboEngineSearchM2M.C
 $(locateModels)/engineSearchIB/engineSearchIB.C
-
+$(locateModels)/hyperEngineSearch/hyperEngineSearch.C
+$(locateModels)/ijkSearch/ijkSearch.C
 
 $(meshMotionModels)/meshMotionModel/meshMotionModel.C
 $(meshMotionModels)/meshMotionModel/newMeshMotionModel.C
 $(meshMotionModels)/noMeshMotion/noMeshMotion.C
+$(meshMotionModels)/DEMdrivenMeshMotion/DEMdrivenMeshMotion.C
 
 $(momCoupleModels)/momCoupleModel/momCoupleModel.C
 $(momCoupleModels)/momCoupleModel/newMomCoupleModel.C
 $(momCoupleModels)/explicitCouple/explicitCouple.C
+$(momCoupleModels)/explicitCoupleSource/explicitCoupleSource.C
 $(momCoupleModels)/implicitCouple/implicitCouple.C
 $(momCoupleModels)/noCouple/noCouple.C
 
-$(regionModels)/regionModel/regionModel.C
-$(regionModels)/regionModel/newRegionModel.C
-$(regionModels)/allRegion/allRegion.C
-
 $(dataExchangeModels)/dataExchangeModel/dataExchangeModel.C
 $(dataExchangeModels)/dataExchangeModel/newDataExchangeModel.C
 $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
@@ -95,11 +120,13 @@ $(dataExchangeModels)/twoWayFiles/twoWayFiles.C
 $(dataExchangeModels)/noDataExchange/noDataExchange.C
 $(dataExchangeModels)/twoWayMPI/twoWayMPI.C
 $(dataExchangeModels)/twoWayM2M/twoWayM2M.C
+$(dataExchangeModels)/twoWayMany2Many/twoWayMany2Many.C
 
 $(averagingModels)/averagingModel/averagingModel.C
 $(averagingModels)/averagingModel/newAveragingModel.C
 $(averagingModels)/dilute/dilute.C
 $(averagingModels)/dense/dense.C
+$(averagingModels)/denseBiDi/denseBiDi.C
 
 $(clockModels)/clockModel/clockModel.C
 $(clockModels)/clockModel/newClockModel.C
@@ -108,6 +135,7 @@ $(clockModels)/noClock/noClock.C
 
 $(liggghtsCommandModels)/liggghtsCommandModel/liggghtsCommandModel.C
 $(liggghtsCommandModels)/liggghtsCommandModel/newLiggghtsCommandModel.C
+$(liggghtsCommandModels)/colorParticles/colorParticles.C
 $(liggghtsCommandModels)/execute/execute.C
 $(liggghtsCommandModels)/runLiggghts/runLiggghts.C
 $(liggghtsCommandModels)/writeLiggghts/writeLiggghts.C
@@ -117,5 +145,6 @@ $(smoothingModels)/smoothingModel/smoothingModel.C
 $(smoothingModels)/smoothingModel/newSmoothingModel.C
 $(smoothingModels)/noSmoothing/noSmoothing.C
 $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C
+$(smoothingModels)/localPSizeDiffSmoothing/localPSizeDiffSmoothing.C
 
-LIB = $(FOAM_USER_LIBBIN)/lib$(CFDEM_LIB_NAME)
+LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME)

src/lagrangian/cfdemParticle/Make/options

@@ -14,6 +14,7 @@ EXE_INC = \
     -I$(LIB_SRC)/OpenFOAM/containers/HashTables/labelHashSet \
     -I$(CFDEM_LIGGGHTS_SRC_DIR) \
     -I$(CFDEM_M2MLIB_PATH) \
+    -I$(CFDEM_Many2ManyLIB_PATH) \
     -I$(CFDEM_SRC_DIR)/cfdTools \
 
 LIB_LIBS = \
@@ -28,4 +29,10 @@ LIB_LIBS = \
     -L$(CFDEM_LIGGGHTS_SRC_DIR) \
     -l$(CFDEM_LIGGGHTS_LIB_NAME) \
     -L$(CFDEM_M2MLIB_PATH) \
-    -lcouple
+    -lcouple \
+    -L$(CFDEM_Many2ManyLIB_PATH) \
+    -lcoupleMany2Many \
+
+/* add -I$(CFDEM_POEMSLIB_PATH) \ to EXE_INC */
+/*    -L$(CFDEM_POEMSLIB_PATH) \ */
+/*    -lpoems   */

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -1,4 +1,4 @@
-word CFDEMversion="cfdem-2.6.3";
+word CFDEMversion="cfdem-2.6.4";
 word compatibleLIGGGHTSversion="3.0.1";
 word OFversion="2.2.x-commit-61b850bc107bdd60bbf1bf9a6417b9faf701d128";
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -73,6 +73,7 @@ Foam::cfdemCloud::cfdemCloud
             IOobject::NO_WRITE
         )
     ),
+    solveFlow_(true),
     verbose_(false),
     ignore_(false),
     modelType_(couplingProperties_.lookup("modelType")),
@@ -95,8 +96,10 @@ Foam::cfdemCloud::cfdemCloud
     momCoupleModels_(couplingProperties_.lookup("momCoupleModels")),
     liggghtsCommandModelList_(liggghtsCommandDict_.lookup("liggghtsCommandModels")),
     turbulenceModelType_(couplingProperties_.lookup("turbulenceModelType")),
+    cg_(1.),
     cgOK_(true),
     impDEMdrag_(false),
+    imExSplitFactor_(1.0),
     useDDTvoidfraction_(false),
     ddtVoidfraction_
     (   
@@ -212,7 +215,10 @@ Foam::cfdemCloud::cfdemCloud
     #include "versionInfo.H"
 
     Info << "If BC are important, please provide volScalarFields -imp/expParticleForces-" << endl;
-
+    if (couplingProperties_.found("solveFlow"))
+        solveFlow_=Switch(couplingProperties_.lookup("solveFlow"));
+    if (couplingProperties_.found("imExSplitFactor"))
+        imExSplitFactor_ = readScalar(couplingProperties_.lookup("imExSplitFactor"));
     if (couplingProperties_.found("verbose")) verbose_=true;
     if (couplingProperties_.found("ignore")) ignore_=true;
     if (turbulenceModelType_=="LESProperties")
@@ -259,7 +265,7 @@ Foam::cfdemCloud::cfdemCloud
     }
 
     dataExchangeM().setCG();
-    if (!cgOK_ && forceM(0).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  * * * * * * * * * * * * * * //
@@ -505,11 +511,47 @@ bool Foam::cfdemCloud::evolve
             if(verbose_) Info << "setvoidFraction done." << endl;
             clockM().stop("setvoidFraction");
 
-            // set particles velocity field
+            // set average particles velocity field
             clockM().start(20,"setVectorAverage");
             setVectorAverages();
-            clockM().stop("setVectorAverage");
 
+
+            //Smoothen "next" fields            
+            smoothingM().dSmoothing();
+            smoothingM().smoothen(voidFractionM().voidFractionNext());
+            smoothingM().smoothenReferenceField(averagingM().UsNext());
+            
+            clockM().stop("setVectorAverage");
+        }
+        
+        //============================================
+        //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;
+        clockM().start(24,"interpolateEulerFields");
+
+        // update voidFractionField
+        alpha = voidFractionM().voidFractionInterp();
+        if(dataExchangeM().couplingStep() < 2)
+        {
+            alpha.oldTime() = alpha; // supress volume src
+            alpha.oldTime().correctBoundaryConditions();
+        }
+        alpha.correctBoundaryConditions();
+
+        // calc ddt(voidfraction)
+        calcDdtVoidfraction(alpha);
+
+        // update mean particle velocity Field
+        Us = averagingM().UsInterp();
+        Us.correctBoundaryConditions();
+
+        clockM().stop("interpolateEulerFields");
+        //============================================
+
+        if(doCouple)
+        {
             // set particles forces
             clockM().start(21,"setForce");
             if(verbose_) Info << "- setForce(forces_)" << endl;
@@ -530,32 +572,7 @@ bool Foam::cfdemCloud::evolve
             giveDEMdata();
             clockM().stop("giveDEMdata");
         }//end dataExchangeM().couple()
-        Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;
 
-        clockM().start(24,"interpolateEulerFields");
-        // update smoothing model
-        smoothingM().dSmoothing();
-
-        //============================================
-        // update voidFractionField V1
-        alpha = voidFractionM().voidFractionInterp();
-        smoothingM().smoothen(alpha);
-        if(dataExchangeM().couplingStep() < 2)
-        {
-            alpha.oldTime() = alpha; // supress volume src
-            alpha.oldTime().correctBoundaryConditions();
-        }
-        alpha.correctBoundaryConditions();
-
-        // calc ddt(voidfraction)
-        calcDdtVoidfraction(alpha);
-
-        // update particle velocity Field
-        Us = averagingM().UsInterp();
-        smoothingM().smoothenReferenceField(Us);
-        Us.correctBoundaryConditions();
-        //============================================
-        clockM().stop("interpolateEulerFields");
 
         if(verbose_){
             #include "debugInfo.H"

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -88,6 +88,8 @@ protected:
 
     IOdictionary liggghtsCommandDict_;
 
+    Switch solveFlow_;
+
     bool verbose_;
 
     bool ignore_;
@@ -132,10 +134,14 @@ protected:
 
     const word turbulenceModelType_;
 
+    mutable scalar cg_;
+
     bool cgOK_;
 
     bool impDEMdrag_;
 
+    mutable scalar imExSplitFactor_;
+
     bool useDDTvoidfraction_;
 
     mutable volScalarField ddtVoidfraction_;
@@ -235,10 +241,20 @@ public:
 
        label liggghtsCommandModelIndex(word);
 
+       inline void setCG(double) const;
+
+       inline const scalar& cg() const;
+
+       inline const bool& impDEMdrag() const;
+
+       inline const scalar& imExSplitFactor() const;
+
        inline const bool& ignore() const;
 
        inline const fvMesh& mesh() const;
 
+       inline bool solveFlow() const;
+
        inline bool verbose() const;
 
        inline const IOdictionary& couplingProperties() const;
@@ -281,6 +297,10 @@ public:
        virtual inline int ** particleTypes() const {return NULL;};
        virtual label particleType(label index) const {return -1;};
 
+       //access to the particles' orientation information
+       virtual inline double ** exArray() const {return NULL;};
+       virtual vector ex(int) const {return Foam::vector(0,0,0);};
+
        inline int numberOfParticles() const;
 
        inline bool numberOfParticlesChanged() const;

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -44,6 +44,27 @@ 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 scalar& cfdemCloud::imExSplitFactor() const
+{
+	    return imExSplitFactor_;
+};
+
+inline const scalar& cfdemCloud::cg() const
+{ 
+    return cg_;
+};
+
 inline const bool& cfdemCloud::ignore() const
 {
     return ignore_;
@@ -54,6 +75,11 @@ inline const fvMesh& cfdemCloud::mesh() const
     return mesh_;
 }
 
+inline bool cfdemCloud::solveFlow() const
+{
+    return bool(solveFlow_);
+}
+
 inline bool cfdemCloud::verbose() const
 {
     return verbose_;

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

@@ -1,5 +1,5 @@
-#define version23  // currently being tested
+//#define version23  // currently being tested
-//#define version22  // currently being used
+#define version22  // currently being used
 //#define version21
 //#define version16ext
 //#define version15
@@ -8,7 +8,7 @@
 //#define comp   // if comp is on - you must use Make/options_comp!
 
 //define multi sphere
-#define multisphere
+//#define multisphere
 
 // features of 2.1 work also in 2.3
 #if defined(version23)

src/lagrangian/cfdemParticle/etc/bashrc

@@ -43,6 +43,10 @@ export CFDEM_LIGGGHTS_LIB_NAME=lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME
 #- CFDEM lib name
 export CFDEM_LIB_NAME=lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
 
+#- LMP Many2Many lib path and makefile
+export CFDEM_Many2ManyLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
+export CFDEM_Many2ManyLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME
+
 #- LMP M2M lib path and makefile
 export CFDEM_M2MLIB_PATH=$CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/library
 export CFDEM_M2MLIB_MAKEFILENAME=$CFDEM_LIGGGHTS_MAKEFILE_NAME
@@ -135,9 +139,6 @@ alias cfdemCompCFDEMsol='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compil
 #- shortcut to compile CFDEMcoupling utilities
 alias cfdemCompCFDEMuti='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh'
 
-#- shortcut to compile couple library
-alias cfdemCompM2M='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileM2Mlib.sh'
-
 #- shortcut to test basic tutorials
 alias cfdemTestTUT='bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/testTutorials.sh'
 
@@ -152,6 +153,10 @@ export -f cfdemLiggghts
 cfdemLiggghtsPar() { mpirun -np $2 $CFDEM_LIGGGHTS_SRC_DIR/lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME < $1; }
 export -f cfdemLiggghtsPar
 
+#- shortcut to open files including a pattern
+cfdemGrep() { grep -rl "$1" ./* | xargs gedit; }
+export -f cfdemGrep
+
 # check if the run directory exists
 if [ -d "$CFDEM_PROJECT_USER_DIR" ]; then
     :

src/lagrangian/cfdemParticle/etc/compileLIGGGHTS_lib.sh

@@ -15,7 +15,6 @@ logDir="log"
 cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc
 mkdir -p $logDir
 
-
 #================================================================================#
 # compile src
 #================================================================================#

src/lagrangian/cfdemParticle/etc/cshrc

@@ -43,6 +43,10 @@ setenv CFDEM_LIGGGHTS_LIB_NAME lmp_$CFDEM_LIGGGHTS_MAKEFILE_NAME
 #- CFDEM lib name
 setenv CFDEM_LIB_NAME lagrangianCFDEM-$CFDEM_VERSION-$WM_PROJECT_VERSION
 
+#- LMP Many2Many lib path and makefile
+setenv CFDEM_Many2ManyLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library
+setenv CFDEM_Many2ManyLIB_MAKEFILENAME $CFDEM_LIGGGHTS_MAKEFILE_NAME
+
 #- LMP M2M lib path
 setenv CFDEM_M2MLIB_PATH $CFDEM_SRC_DIR/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayM2M/library
 setenv CFDEM_M2MLIB_MAKEFILENAME $CFDEM_LIGGGHTS_MAKEFILE_NAME
@@ -133,9 +137,6 @@ alias cfdemCompCFDEMsol 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compil
 #- shortcut to compile CFDEMcoupling utilities
 alias cfdemCompCFDEMuti 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileCFDEMcoupling_uti.sh'
 
-#- shortcut to compile couple library
-alias cfdemCompM2M 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/compileM2Mlib.sh'
-
 #- shortcut to test basic tutorials
 alias cfdemTestTUT 'bash $CFDEM_SRC_DIR/lagrangian/cfdemParticle/etc/testTutorials.sh'
 

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

@@ -1 +1,11 @@
 lagrangian/cfdemParticle/dir
+
+#====================================================='
+#- RADL
+fvOptions/dir
+cylPorousMedia/dir
+
+#====================================================='
+#- other
+finiteVolume/dir
+

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

@@ -1,4 +1,13 @@
-cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
+cfdemSolverPimpleImEx/dir
+cfdemSolverIBInterLubrication/dir
+cfdemSolverIBScalar/dir
+cfdemSolverInterDyM/dir
+cfdemSolverInterDyMPC/dir
+cfdemSolverBubble/dir
+cfdemSolverPisoMS/dir
+cfdemSolverPimpleDyM_22x/dir
+cfdemSolverPimpleDyMMS_22x/dir
+cfdemSolverPimpleDyMScalar_22x/dir

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

@@ -7,15 +7,51 @@
 #===================================================================#
 
 cfdemSolverPiso/settlingTestMPI/dir
-
 cfdemSolverPiso/ErgunTestMPI/dir
-
 cfdemSolverPiso/ErgunTestMPI_cgs/dir
-
 cfdemSolverPiso/ErgunTestMPI_restart/dir
-
 cfdemSolverIB/twoSpheresGlowinskiMPI/dir
-
 cfdemSolverPisoScalar/packedBedTemp/dir
 
-cfdemSolverPiso/ErgunTestCG/dir
+#===================================================================#
+# RADL
+cfdemSolverPimpleImEx/settlingTestMPI/dir
+cfdemSolverPimpleImEx/ErgunTestMPI/dir
+#cfdemSolverPimpleImEx/crossFlow/dir
+#cfdemSolverIB/periodicCase/dir
+#cfdemSolverIB/cfdemIBPeriodicCubicalBox_fullyPeriodic/dir
+#cfdemSolverIBInterLubrication/twoCoatedParticlesRelMotion_smallTest/dir
+#cfdemSolverIBScalar/cfdemIBPeriodicCubicalBoxScalar/dir
+
+#===================================================================#
+# NesteJacobs
+#Projects/Neste/cfdemSolverBubble/3pFBreactor/dir
+#Projects/Neste/cfdemSolverInterDyM/3pFBreactor/dir
+
+#===================================================================#
+# not in release:
+
+#cfdemSolverPiso/settlingTestBigParticleMPI/dir
+#cfdemSolverPiso/ErgunTestCG/dir
+#cfdemSolverPiso/ErgunTestM2M/dir
+#cfdemSolverPiso/HopperEmptying/dir
+
+cfdemSolverPimpleDyM/ErgunTestMPI/dir
+
+#cfdemSolverPisoMS/settlingTestMPI/dir
+#cfdemSolverPisoMS/ErgunTestMPI/dir
+
+#cfdemSolverInterDyM/twoPhaseSettlingTest/dir
+#cfdemSolverInterDyM/ErgunTestMPI/dir
+#cfdemSolverInterDyM/granularPiston/dir
+#cfdemSolverInterDyM/sugarNcoffee/dir
+
+#cfdemSolverBubble/ErgunTestMPI_pureLiquid/dir
+
+#- these examples are already designed for 2.3.x
+#cfdemSolverInterDyMPC/sugarNcoffee/dir
+#cfdemSolverInterDyMPC/granularPiston/dir
+#cfdemSolverInterDyMPC/meltingPot/dir
+
+
+

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

@@ -45,8 +45,10 @@ defineRunTimeSelectionTable(IOModel, dictionary);
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void IOModel::dumpDEMdata() const
+int IOModel::dumpDEMdata() const
-{}
+{
+    return -1;
+}
 
 fileName IOModel::createTimeDir(fileName path) const
 {
@@ -64,6 +66,60 @@ fileName IOModel::createLagrangianDir(fileName path) const
     return cfdemCloudDirPath;
 }
 
+fileName IOModel::buildFilePath(word dirName) const
+{
+    // create file structure
+	fileName path("");
+    if(parOutput_)
+    {
+    	path=fileName(particleCloud_.mesh().time().path()/particleCloud_.mesh().time().timeName()/dirName/"particleCloud");
+    	mkDir(path,0777);
+    } else
+    {
+		path=fileName("."/dirName);
+    	mkDir(path,0777);
+    	mkDir(fileName(path/"constant"),0777);
+    	OFstream* stubFile = new OFstream(fileName(path/"particles.foam"));
+    	delete stubFile;
+        }
+    return path;
+}
+
+void IOModel::streamDataToPath(fileName path, double** array,int nPProc,word name,word type,word className,word finaliser) const
+{
+    vector vec;
+    OFstream* fileStream = new OFstream(fileName(path/name));
+    *fileStream << "FoamFile\n";
+    *fileStream << "{version 2.0; format ascii;class "<< className << "; location 0;object  "<< name <<";}\n";
+    *fileStream << nPProc <<"\n";
+
+
+    if(type!="origProcId")*fileStream << "(\n";
+    else if(type=="origProcId")
+    {
+        if(nPProc>0) *fileStream <<"{0}"<< "\n";
+        else *fileStream <<"{}"<< "\n";
+    }
+
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+        if (particleCloud_.cellIDs()[index][0] > -1) // particle Found
+        {
+            if (type=="scalar"){
+                *fileStream << array[index][0] << " \n";
+            }else if (type=="position" || type=="vector"){
+                for(int i=0;i<3;i++) vec[i] = array[index][i];
+                *fileStream <<"( "<< vec[0] <<" "<<vec[1]<<" "<<vec[2]<<" ) "<< finaliser << " \n";
+            }else if (type=="label"){
+                *fileStream << index << finaliser << " \n";
+            }
+        }
+    }
+
+    if(type!="origProcId")*fileStream << ")\n";
+    delete fileStream;
+}
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
@@ -75,8 +131,18 @@ IOModel::IOModel
 :
     dict_(dict),
     particleCloud_(sm),
-    time_(sm.mesh().time())
+    time_(sm.mesh().time()),
-{}
+    parOutput_(true)
+{
+	if (
+            particleCloud_.dataExchangeM().myType()=="oneWayVTK" ||
+            dict_.found("serialOutput")
+       )
+    {
+        parOutput_=false;
+        Warning << "IO model is in serial write mode, only data on proc 0 is written" << endl;
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

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

@@ -63,6 +63,8 @@ protected:
 
         const Time& time_;
 
+        bool parOutput_;
+
 public:
 
     //- Runtime type information
@@ -109,11 +111,15 @@ public:
 
     // Member Functions
 
-    virtual void dumpDEMdata() const;
+    virtual int dumpDEMdata() const;
 
     fileName createTimeDir(fileName) const;
 
     fileName createLagrangianDir(fileName) const;
+
+    fileName buildFilePath(word) const;
+
+    void streamDataToPath(fileName,double**,int,word,word type,word className,word finaliser) const;
 };
 
 

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

@@ -64,19 +64,9 @@ basicIO::basicIO
     //propsDict_(dict.subDict(typeName + "Props")),
     dirName_("lagrangian"),
     path_("dev/null"),
-    parOutput_(true),
     nPProc_(-1),
     lagPath_("dev/null")
 {
-	if (
-            particleCloud_.dataExchangeM().myType()=="oneWayVTK" ||
-            dict_.found("serialOutput")
-       )
-    {
-        parOutput_=false;
-        Warning << "IO model is in serial write mode, only data on proc 0 is written" << endl;
-    }
-
     //if (propsDict_.found("dirName")) dirName_=word(propsDict_.lookup("dirName"));
     path_ = buildFilePath(dirName_);
 }
@@ -92,7 +82,7 @@ basicIO::~basicIO()
 
 // Public Member Functions
 
-void basicIO::dumpDEMdata() const
+int basicIO::dumpDEMdata() const
 {
     if (time_.outputTime())
     {
@@ -112,67 +102,16 @@ void basicIO::dumpDEMdata() const
         nPProc_=count;
         
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.positions(), "positions","vector","Cloud<passiveParticle>","0");
+        streamDataToPath(lagPath_, particleCloud_.positions(),nPProc_,"positions","vector","Cloud<passiveParticle>","0");
-        streamDataToPath(lagPath_, particleCloud_.velocities(), "v","vector","vectorField","");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),nPProc_,"v","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.radii(), "r","scalar","scalarField","");
+        streamDataToPath(lagPath_, particleCloud_.radii(),nPProc_,"r","scalar","scalarField","");
     }
+    return nPProc_;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Private Member Functions
 
-fileName basicIO::buildFilePath(word dirName) const
-{
-    // create file structure
-	fileName path("");
-    if(parOutput_)
-    {
-    	path=fileName(particleCloud_.mesh().time().path()/particleCloud_.mesh().time().timeName()/dirName/"particleCloud");
-    	mkDir(path,0777);
-    } else
-    {
-		path=fileName("."/dirName);
-    	mkDir(path,0777);
-    	mkDir(fileName(path/"constant"),0777);
-    	OFstream* stubFile = new OFstream(fileName(path/"particles.foam"));
-    	delete stubFile;
-        }
-    return path;
-}
 
-void basicIO::streamDataToPath(fileName path, double** array,word name,word type,word className,word finaliser) const
-{
-    vector vec;
-    OFstream* fileStream = new OFstream(fileName(path/name));
-    *fileStream << "FoamFile\n";
-    *fileStream << "{version 2.0; format ascii;class "<< className << "; location 0;object  "<< name <<";}\n";
-    *fileStream << nPProc_ <<"\n";
-    //*fileStream << "(\n";
-
-    if(type!="origProcId")*fileStream << "(\n";
-    else if(type=="origProcId")
-    {
-        if(nPProc_>0) *fileStream <<"{0}"<< "\n";
-        else *fileStream <<"{}"<< "\n";
-    }
-
-    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
-    {
-        if (particleCloud_.cellIDs()[index][0] > -1) // particle Found
-        {
-            if (type=="scalar"){
-                *fileStream << array[index][0] << " \n";
-            }else if (type=="position" || type=="vector"){
-                for(int i=0;i<3;i++) vec[i] = array[index][i];
-                *fileStream <<"( "<< vec[0] <<" "<<vec[1]<<" "<<vec[2]<<" ) "<< finaliser << " \n";
-            }else if (type=="label"){
-                *fileStream << index << finaliser << " \n";
-            }
-        }
-    }
-    //*fileStream << ")\n";
-    if(type!="origProcId")*fileStream << ")\n";
-    delete fileStream;
-}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -64,18 +64,13 @@ private:
 
     fileName path_;
 
-    bool parOutput_;
-
     mutable int nPProc_;
 
     // Member Functions
-    fileName buildFilePath(word) const;
 
 protected:
     mutable fileName lagPath_;
 
-    void streamDataToPath(fileName,double**,word,word type,word className,word finaliser) const;
-
 public:
 
     //- Runtime type information
@@ -97,7 +92,7 @@ public:
 
 
     // Member Functions
-    void dumpDEMdata() const;
+    int dumpDEMdata() const;
 };
 
 

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

@@ -75,18 +75,21 @@ sophIO::~sophIO()
 
 // Public Member Functions
 
-void sophIO::dumpDEMdata() const
+int sophIO::dumpDEMdata() const
 {
+    int npProcs(-1);
+
     if (time_.outputTime())
     {
-        basicIO::dumpDEMdata();
+        npProcs=basicIO::dumpDEMdata();
 
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.voidfractions(), "voidfractions","scalar","scalarField","");
+        streamDataToPath(lagPath_, particleCloud_.voidfractions(),npProcs,"voidfractions","scalar","scalarField","");
-        streamDataToPath(lagPath_, particleCloud_.impForces(), "impForces","vector","vectorField","");
+        streamDataToPath(lagPath_, particleCloud_.impForces(),npProcs,"impForces","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.expForces(), "expForces","vector","vectorField","");
+        streamDataToPath(lagPath_, particleCloud_.expForces(),npProcs,"expForces","vector","vectorField","");
-        streamDataToPath(lagPath_, particleCloud_.DEMForces(), "DEMForces","vector","vectorField","");
+        streamDataToPath(lagPath_, particleCloud_.DEMForces(),npProcs,"DEMForces","vector","vectorField","");
     }
+    return npProcs;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Private Member Functions

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

@@ -81,7 +81,7 @@ public:
 
 
     // Member Functions
-    void dumpDEMdata() const;
+    int dumpDEMdata() const;
 };
 
 

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

@@ -73,16 +73,20 @@ trackIO::~trackIO()
 
 // Public Member Functions
 
-void trackIO::dumpDEMdata() const
+int trackIO::dumpDEMdata() const
 {
+    int npProcs(-1);
+
     if (time_.outputTime())
     {
-        sophIO::dumpDEMdata();
+        npProcs = sophIO::dumpDEMdata();
 
         // stream data to file
-        streamDataToPath(lagPath_, particleCloud_.velocities(), "origId","label","labelField","");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origId","label","labelField","");
-        streamDataToPath(lagPath_, particleCloud_.velocities(), "origProcId","origProcId","labelField","");
+        streamDataToPath(lagPath_, particleCloud_.velocities(),npProcs,"origProcId","origProcId","labelField","");
     }
+
+    return npProcs;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Private Member Functions

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

@@ -80,7 +80,7 @@ public:
 
 
     // Member Functions
-    void dumpDEMdata() const;
+    int dumpDEMdata() const;
 };
 
 

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

@@ -245,7 +245,7 @@ void averagingModel::setDSauter
     scalar radiusPow3(-1);
     scalar volume(-1);
 
-    scalar scale_ = particleCloud_.forceM(0).cg(); //scaling of parcel vs. primary particle diameter
+    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)
     volScalarField riPower2
     (

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

@@ -128,7 +128,6 @@ twoWayM2M::twoWayM2M
     lmp2foam_vec_ = NULL;
     foam2lmp_vec_ = NULL;
     foam2lmp_ = NULL;
-    foam2lmp_ = NULL;
     nlocal_lammps_ = -1;
     id_lammps_ = NULL;
     id_lammpsVec_ = NULL;
@@ -170,7 +169,6 @@ twoWayM2M::~twoWayM2M()
     delete lmp2foam_vec_;
     delete foam2lmp_vec_;
     delete foam2lmp_;
-    delete foam2lmp_;
     delete lmp;
 }
 
@@ -264,6 +262,7 @@ void twoWayM2M::giveData
         //==================
     }else if( type == "scalar-atom" )
     {
+        Warning << "LIGGGHTS not ready for use of impleDEM and Many2Many" << endl;
         double *tmp_=NULL;
         LAMMPS_NS::Fix *fix = NULL;
         fix = lmp->modify->find_fix_property(charName,"property/atom","scalar",0,0,"cfd coupling",false);
@@ -320,7 +319,7 @@ void Foam::twoWayM2M::allocateArray
         for (int j = 0; j < width; j++)
             array[i][j] = initVal;
 }
-void Foam::twoWayM2M::destroy(double** array,int len) const
+void inline Foam::twoWayM2M::destroy(double** array,int len) const
 {
     // destroy array[i] first?
     lmp->memory->destroy(array);
@@ -356,7 +355,7 @@ void Foam::twoWayM2M::allocateArray
         for (int j = 0; j < width; j++)
             array[i][j] = initVal;
 }
-void Foam::twoWayM2M::destroy(int** array,int len) const
+void inline Foam::twoWayM2M::destroy(int** array,int len) const
 {
     // destroy array[i] first?
     lmp->memory->destroy(array);
@@ -370,7 +369,7 @@ void Foam::twoWayM2M::allocateArray(double*& array, double initVal, int length)
     for (int i = 0; i < len; i++)
         array[i] = initVal;
 }
-void Foam::twoWayM2M::destroy(double* array) const
+void inline Foam::twoWayM2M::destroy(double* array) const
 {
     lmp->memory->destroy(array);
 }
@@ -383,7 +382,7 @@ void Foam::twoWayM2M::allocateArray(int*& array, int initVal, int length) const
     for (int i = 0; i < len; i++)
         array[i] = initVal;
 }
-void Foam::twoWayM2M::destroy(int* array) const
+void inline Foam::twoWayM2M::destroy(int* array) const
 {
     lmp->memory->destroy(array);
 }
@@ -409,7 +408,7 @@ bool Foam::twoWayM2M::couple() const
 
                 if(particleCloud_.liggghtsCommand()[i]().runCommand(couplingStep()))
                 {
-                    const char* command = particleCloud_.liggghtsCommand()[i]().command();
+                    const char* command = particleCloud_.liggghtsCommand()[i]().command(0);
                     Info << "Executing command: '"<< command <<"'"<< endl;
                     lmp->input->one(command);
                 }
@@ -485,6 +484,8 @@ void Foam::twoWayM2M::syncIDs() const
 
     int*  id_lammpsSync=NULL;
     double** pos_lammpsSync=NULL;
+	bool pos_lammps_alloc_flag=false;
+	bool id_lammps_alloc_flag=false;
 
     if(firstRun_ || particleLost_)
     {
@@ -503,9 +504,11 @@ void Foam::twoWayM2M::syncIDs() const
             for (int j=0;j<3;j++)
                 id_lammpsVec_[i*3+j] = id_lammps_[i]*3+j;
         }
-
+		destroy(pos_lammps_,0);
-        // get access to "x"
+		pos_lammps_=NULL;
         pos_lammps_ = (double **) lammps_extract_atom(lmp,"x");
+		pos_lammps_alloc_flag=false;
+		id_lammps_alloc_flag=false;
     }
     else
     {
@@ -532,6 +535,8 @@ void Foam::twoWayM2M::syncIDs() const
         // map data according to last TS
         id_lammps_=NULL;
         allocateArray(id_lammps_,-1.,nlocal_foam_);
+		id_lammps_alloc_flag=true;
+
         allocateArray(tmpI_,-1.,nlocal_foam_);
         lmp2foam_->exchange(id_lammpsSync, tmpI_);
         for(int i=0;i<nlocal_foam_;i++)
@@ -544,6 +549,7 @@ void Foam::twoWayM2M::syncIDs() const
         lmp2foam_vec_->exchange(pos_lammpsSync ? pos_lammpsSync[0] : NULL, tmp_);
 
         allocateArray(pos_lammps_,0,3,nlocal_foam_);
+		pos_lammps_alloc_flag=true;
         for(int i=0;i<nlocal_foam_;i++)
             for(int j=0;j<3;j++)
                 pos_lammps_[i][j]=tmp_[i*3+j];
@@ -554,7 +560,8 @@ void Foam::twoWayM2M::syncIDs() const
     particleCloud_.clockM().stop("recv_DEM_ids");
 
     particleCloud_.clockM().start(6,"locateParticle()");
-    locateParticle(id_lammpsSync);
+    locateParticle(id_lammpsSync, id_lammps_alloc_flag);
+	id_lammps_alloc_flag=true;
     particleCloud_.clockM().stop("locateParticle()");
 
 //MPI_Barrier(MPI_COMM_WORLD);
@@ -619,8 +626,11 @@ void Foam::twoWayM2M::syncIDs() const
         foam2lmp_vec_->setup(nlocal_foam_*3,id_foamVec_,nlocal_lammps_*3,id_lammpsVec_);
         foam2lmp_->setup(nlocal_foam_,id_foam_,nlocal_lammps_,id_lammpsSync);
     }
+    if (id_lammps_alloc_flag) destroy(id_lammps_);
     id_lammps_=NULL;    // free pointer from LIG
+    destroy(id_lammpsSync);
     id_lammpsSync=NULL; // free pointer from LIG
+    if (pos_lammps_alloc_flag) destroy(pos_lammps_,0);
     pos_lammps_ = NULL; // free pointer from LIG
 
     particleCloud_.clockM().stop("setup_Comm");
@@ -629,7 +639,7 @@ void Foam::twoWayM2M::syncIDs() const
 //Info << "update communication - done." << endl;
 }
 
-void Foam::twoWayM2M::locateParticle(int* id_lammpsSync) const
+void Foam::twoWayM2M::locateParticle(int* id_lammpsSync, bool id_lammps_alloc_flag) const
 {
 #if defined(version21)
 
@@ -918,6 +928,9 @@ void Foam::twoWayM2M::locateParticle(int* id_lammpsSync) const
             }
             i++;
         }
+		Foam::dataExchangeModel::destroy(id_foam_nowhere_all);
+		id_foam_nowhere_all=NULL;
+		if (id_lammps_alloc_flag) destroy(id_lammps_);
 
         // make cpy
         id_lammps_=NULL;

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

@@ -188,22 +188,22 @@ public:
         // double **
         void allocateArray(double**&, double, int, int) const;
         void allocateArray(double**&, double, int,const char* ="nparticles") const;
-        void destroy(double**,int) const;
+        void inline destroy(double**,int) const;
         //============
         // int **
         void allocateArray(int**&, int, int, int) const;
         void allocateArray(int**&, int, int,const char* ="nparticles") const;
-        void destroy(int**,int) const;
+        void inline destroy(int**,int) const;
         //==============
 
         //==============
         // double *
         void allocateArray(double*&, double, int) const;
-        void destroy(double*) const;
+        void inline destroy(double*) const;
         //==============
         // int *
         void allocateArray(int*&, int, int) const;
-        void destroy(int*) const;
+        void inline destroy(int*) const;
         //==============
 
         bool couple() const;
@@ -212,9 +212,9 @@ public:
         int getNumberOfClumps() const;
 
         void syncIDs() const;
-        void locateParticle(int*) const;
+        void locateParticle(int*, bool) const;
         word myType() const{return typeName; };
-        void setCG() const { particleCloud_.forceM(0).setCG(lmp->force->cg()); };
+        void setCG() const { particleCloud_.setCG(lmp->force->cg()); };
 };
 
 

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

@@ -327,23 +327,25 @@ bool Foam::twoWayMPI::couple() const
             }
 
             // models with exact timing exists
+            label commandLines(0);
             if(exactTiming)
             {
                 // extension for more liggghtsCommands active the same time:
                 //    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];          
                     particleCloud_.liggghtsCommand()[runComNr]().set(DEMstepsToInterrupt[j]);
-                    const char* command = particleCloud_.liggghtsCommand()[runComNr]().command();
+                    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();
+                    command = particleCloud_.liggghtsCommand()[lcModel[j]]().command(0);
                     Info << "Executing command: '"<< command <<"'"<< endl;
                     lmp->input->one(command);
                 }
@@ -352,7 +354,7 @@ bool Foam::twoWayMPI::couple() const
                 if(particleCloud_.liggghtsCommand()[runComNr]().runCommand(couplingStep()))
                 {
                     particleCloud_.liggghtsCommand()[runComNr]().set(couplingInterval() - DEMstepsRun);
-                    const char* command = particleCloud_.liggghtsCommand()[runComNr]().command();
+                    const char* command = particleCloud_.liggghtsCommand()[runComNr]().command(0);
                     Info << "Executing run command: '"<< command <<"'"<< endl;
                     lmp->input->one(command);
                 }
@@ -366,9 +368,13 @@ bool Foam::twoWayMPI::couple() const
                         particleCloud_.liggghtsCommand()[i]().runCommand(couplingStep())
                     )
                     {
-                        const char* command = particleCloud_.liggghtsCommand()[i]().command();
+                        commandLines=particleCloud_.liggghtsCommand()[i]().commandLines();
-                        Info << "Executing command: '"<< command <<"'"<< endl;
+                        for(int j=0;j<commandLines;j++)
-                        lmp->input->one(command);
+                        {
+                            const char* command = particleCloud_.liggghtsCommand()[i]().command(j);
+                            Info << "Executing command: '"<< command <<"'"<< endl;
+                            lmp->input->one(command);
+                        }
                     }
                 }
             }
@@ -379,9 +385,13 @@ bool Foam::twoWayMPI::couple() const
                 {
                     if(particleCloud_.liggghtsCommand()[i]().runCommand(couplingStep()))
                     {
-                        const char* command = particleCloud_.liggghtsCommand()[i]().command();
+                        commandLines=particleCloud_.liggghtsCommand()[i]().commandLines();
-                        Info << "Executing command: '"<< command <<"'"<< endl;
+                        for(int j=0;j<commandLines;j++)
-                        lmp->input->one(command);
+                        {
+                            const char* command = particleCloud_.liggghtsCommand()[i]().command(j);
+                            Info << "Executing command: '"<< command <<"'"<< endl;
+                            lmp->input->one(command);
+                        }
                     }
                 }
             }

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

@@ -167,7 +167,7 @@ public:
 
         word myType() const{return typeName; };
 
-        void setCG() const { particleCloud_.forceM(0).setCG(lmp->force->cg()); };
+        void setCG() const { particleCloud_.setCG(lmp->force->cg()); };
 };
 
 

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

@@ -71,7 +71,12 @@ DiFeliceDrag::DiFeliceDrag
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
-    interpolation_(false)
+    interpolation_(false),
+    splitImplicitExplicit_(false),
+    UsFieldName_(propsDict_.lookup("granVelFieldName")),
+    UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
+    scaleDia_(1.),
+    scaleDrag_(1.)
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, "diFeliceDrag.logDat");
@@ -89,7 +94,19 @@ DiFeliceDrag::DiFeliceDrag
         Info << "using interpolated value of U." << endl;
         interpolation_=true;
     }
-    particleCloud_.checkCG(false);
+    if (propsDict_.found("splitImplicitExplicit"))
+    {
+        Info << "will split implicit / explicit force contributions." << endl;
+        splitImplicitExplicit_ = true;
+        if(!interpolation_) 
+            Info << "WARNING: will only consider fluctuating particle velocity in implicit / explicit force split!" << endl;
+    }
+    
+    particleCloud_.checkCG(true);
+    if (propsDict_.found("scale"))
+        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+    if (propsDict_.found("scaleDrag"))
+        scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
 }
 
 
@@ -103,6 +120,13 @@ DiFeliceDrag::~DiFeliceDrag()
 
 void DiFeliceDrag::setForce() const
 {
+    if (scaleDia_ > 1)
+        Info << "DiFeliceDrag using scale = " << scaleDia_ << endl;
+    else if (particleCloud_.cg() > 1){
+        scaleDia_=particleCloud_.cg();
+        Info << "DiFeliceDrag using scale from liggghts cg = " << scaleDia_ << endl;
+    }
+    
     // get viscosity field
     #ifdef comp
         const volScalarField nufField = particleCloud_.turbulence().mu() / rho_;
@@ -123,6 +147,11 @@ void DiFeliceDrag::setForce() const
     scalar magUr(0);
     scalar Rep(0);
 	scalar Cd(0);
+	
+	vector UfluidFluct(0,0,0);
+    vector UsFluct(0,0,0);
+    vector dragExplicit(0,0,0);
+  	scalar dragCoefficient(0);
 
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
@@ -136,7 +165,7 @@ void DiFeliceDrag::setForce() const
 
             cellI = particleCloud_.cellIDs()[index][0];
             drag = vector(0,0,0);
-
+            
             if (cellI > -1) // particle Found
             {
                 if(interpolation_)
@@ -158,11 +187,13 @@ void DiFeliceDrag::setForce() const
                 magUr = mag(Ur);
                 Rep = 0;
                 Cd = 0;
+                dragCoefficient = 0;
 
                 if (magUr > 0)
                 {
+                   
                     // calc particle Re Nr
-                    Rep = ds*voidfraction*magUr/(nuf+SMALL);
+                    Rep = ds/scaleDia_*voidfraction*magUr/(nuf+SMALL);
 
                     // calc fluid drag Coeff
                     Cd = sqr(0.63 + 4.8/sqrt(Rep));
@@ -171,24 +202,45 @@ void DiFeliceDrag::setForce() const
                     scalar Xi = 3.7 - 0.65 * exp(-sqr(1.5-log10(Rep))/2);
 
                     // calc particle's drag
-                    drag = 0.125*Cd*rho*M_PI*ds*ds*pow(voidfraction,(2-Xi))*magUr*Ur;
+                    dragCoefficient = 0.125*Cd*rho
-
+                                     *M_PI
+                                     *ds*ds     
+                                     *scaleDia_ 
+                                     *pow(voidfraction,(2-Xi))*magUr
+                                     *scaleDrag_;
                     if (modelType_=="B")
-                        drag /= voidfraction;
+                        dragCoefficient /= voidfraction;
+                        
+                    drag = dragCoefficient*Ur; //total drag force!
+                        
+                    //Split forces
+                    if(splitImplicitExplicit_)
+                    {
+                        UfluidFluct  = Ufluid - U_[cellI];
+                        UsFluct      = Us     - UsField_[cellI];
+                        dragExplicit = dragCoefficient*(UfluidFluct - UsFluct); //explicit part of force
+                    }
+
                 }
 
-                if(verbose_ && index >100 && index <102)
+                if(verbose_ && index >-1 && index <102)
                 {
                     Pout << "index = " << index << endl;
                     Pout << "Us = " << Us << endl;
                     Pout << "Ur = " << Ur << endl;
-                    Pout << "ds = " << ds << endl;
+                    Pout << "ds/scale = " << ds/scaleDia_ << endl;
                     Pout << "rho = " << rho << endl;
                     Pout << "nuf = " << nuf << endl;
                     Pout << "voidfraction = " << voidfraction << endl;
                     Pout << "Rep = " << Rep << endl;
                     Pout << "Cd = " << Cd << endl;
-                    Pout << "drag = " << drag << endl;
+                    Pout << "drag (total) = " << drag << endl;
+                    if(splitImplicitExplicit_)
+                    {
+                        Pout << "UfluidFluct = " << UfluidFluct << endl;
+                        Pout << "UsFluct = " << UsFluct << endl;
+                        Pout << "dragExplicit = " << dragExplicit << endl;
+                    }
                 }
 
                 //Set value fields and write the probe
@@ -205,7 +257,15 @@ void DiFeliceDrag::setForce() const
             }
             // set force on particle
             if(treatExplicit_) for(int j=0;j<3;j++) expForces()[index][j] += drag[j];
-            else  for(int j=0;j<3;j++) impForces()[index][j] += drag[j];
+            else   //implicit treatment, taking explicit force contribution into account
+            {
+               for(int j=0;j<3;j++) 
+               { 
+                    impForces()[index][j] += drag[j] - dragExplicit[j]; //only consider implicit part!
+                    expForces()[index][j] += dragExplicit[j];
+               }
+            }
+            
             for(int j=0;j<3;j++) DEMForces()[index][j] += drag[j];
         }
     //}

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

@@ -76,6 +76,16 @@ private:
 
     bool interpolation_; // use interpolated U field values
 
+    bool splitImplicitExplicit_; // use splitting of implicit and explict force contribution
+
+    word UsFieldName_;   
+
+    const volVectorField& UsField_;  // the average particle velocity field (for implicit/expliti force split)
+
+    mutable scalar scaleDia_;
+
+    mutable scalar scaleDrag_;
+
 public:
 
     //- Runtime type information

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

@@ -71,21 +71,34 @@ GidaspowDrag::GidaspowDrag
     voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     phi_(readScalar(propsDict_.lookup("phi"))),
-    interpolation_(false)
+    interpolation_(false),
+    scaleDia_(1.),
+    scaleDrag_(1.)
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, "gidaspowDrag.logDat");
     particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must  be the force
     particleCloud_.probeM().vectorFields_.append("Urel");
     particleCloud_.probeM().scalarFields_.append("Rep");
-    particleCloud_.probeM().scalarFields_.append("beta");
+    particleCloud_.probeM().scalarFields_.append("betaP");
     particleCloud_.probeM().scalarFields_.append("voidfraction");
     particleCloud_.probeM().writeHeader();
 
     if (propsDict_.found("verbose")) verbose_=true;
     if (propsDict_.found("treatExplicit")) treatExplicit_=true;
     if (propsDict_.found("interpolation")) interpolation_=true;
-    particleCloud_.checkCG(false);
+    if (propsDict_.found("implDEM"))
+    {
+        treatExplicit_=false;
+        implDEM_=true;
+        setImpDEMdrag();
+        Info << "Using implicit DEM drag formulation." << endl;
+    }
+    particleCloud_.checkCG(true);
+    if (propsDict_.found("scale"))
+        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+    if (propsDict_.found("scaleDrag"))
+        scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
 
     Info << "Gidaspow - interpolation switch: " << interpolation_ << endl;
 }
@@ -101,6 +114,13 @@ GidaspowDrag::~GidaspowDrag()
 
 void GidaspowDrag::setForce() const
 {
+    if (scaleDia_ > 1)
+        Info << "Gidaspow using scale = " << scaleDia_ << endl;
+    else if (particleCloud_.cg() > 1){
+        scaleDia_=particleCloud_.cg();
+        Info << "Gidaspow using scale from liggghts cg = " << scaleDia_ << endl;
+    }
+
     // get viscosity field
     #ifdef comp
         const volScalarField nufField = particleCloud_.turbulence().mu() / rho_;
@@ -111,6 +131,8 @@ void GidaspowDrag::setForce() const
     vector position(0,0,0);
     scalar voidfraction(1);
     vector Ufluid(0,0,0);
+    vector drag(0,0,0);
+    label cellI=0;
 
     vector Us(0,0,0);
     vector Ur(0,0,0);
@@ -122,9 +144,9 @@ void GidaspowDrag::setForce() const
     scalar Vs(0);
     scalar localPhiP(0);
 
-    scalar CdMagUrLag(0);    //Cd of the very particle
+    scalar CdMagUrLag(0);       //Cd of the very particle
-    scalar KslLag(0);              //momentum exchange of the very particle (per unit volume)
+    scalar KslLag(0);           //momentum exchange of the very particle (per unit volume)
-    scalar beta(0);              //momentum exchange of the very particle
+    scalar betaP(0);             //momentum exchange of the very particle
 
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
@@ -135,8 +157,12 @@ void GidaspowDrag::setForce() const
     {
         //if(mask[index][0])
         //{
-            vector drag(0,0,0);
+            cellI = particleCloud_.cellIDs()[index][0];
-            label cellI = particleCloud_.cellIDs()[index][0];
+            drag = vector(0,0,0);
+            betaP = 0;
+            Vs = 0;
+            Ufluid =vector(0,0,0);
+            voidfraction=0;
 
             if (cellI > -1) // particle Found
             {
@@ -148,8 +174,8 @@ void GidaspowDrag::setForce() const
                     Ufluid       = UInterpolator_.interpolate(position,cellI);
                     //Ensure interpolated void fraction to be meaningful
                     // Info << " --> voidfraction: " << voidfraction << endl;
-                    if(voidfraction>1.00) voidfraction = 1.0f;
+                    if(voidfraction>1.00) voidfraction = 1.0;
-                    if(voidfraction<0.10) voidfraction = 0.10f;
+                    if(voidfraction<0.10) voidfraction = 0.10;
                 }
                 else
                 {
@@ -171,30 +197,30 @@ void GidaspowDrag::setForce() const
                 //Compute specific drag coefficient (i.e., Force per unit slip velocity and per m³ SUSPENSION)
                 if(voidfraction > 0.8) //dilute
                 {
-                    Rep=ds*voidfraction*magUr/nuf;
+                    Rep=ds/scaleDia_*voidfraction*magUr/nuf;
-                    CdMagUrLag = (24.0*nuf/(ds*voidfraction)) //1/magUr missing here, but compensated in expression for KslLag!
+                    CdMagUrLag = (24.0*nuf/(ds/scaleDia_*voidfraction)) //1/magUr missing here, but compensated in expression for KslLag!
                                  *(scalar(1)+0.15*Foam::pow(Rep, 0.687));
 
                     KslLag = 0.75*(
                                             rho*localPhiP*voidfraction*CdMagUrLag
                                           /
-                                            (ds*Foam::pow(voidfraction,2.65))
+                                            (ds/scaleDia_*Foam::pow(voidfraction,2.65))
                                           );
                 }
                 else  //dense
                 {
                     KslLag = (150*Foam::pow(localPhiP,2)*nuf*rho)/
-                             (voidfraction*ds*ds+SMALL)
+                             (voidfraction*ds/scaleDia_*ds/scaleDia_+SMALL)
                             +
                              (1.75*(localPhiP) * magUr * rho)/
-                             ((ds));
+                             ((ds/scaleDia_));
                 }
 
                 // calc particle's drag coefficient (i.e., Force per unit slip velocity and per m³ PARTICLE)
-                beta = KslLag / localPhiP;
+                betaP = KslLag / localPhiP;
 
                 // calc particle's drag
-                drag = Vs * beta * Ur;
+                drag = Vs * betaP * Ur * scaleDrag_;
 
                 if (modelType_=="B")
                     drag /= voidfraction;
@@ -206,11 +232,13 @@ void GidaspowDrag::setForce() const
                     Pout << "Us = " << Us << endl;
                     Pout << "Ur = " << Ur << endl;
                     Pout << "ds = " << ds << endl;
+                    Pout << "ds/scale = " << ds/scaleDia_ << endl;
+                    Pout << "phi = " << phi_ << endl;
                     Pout << "rho = " << rho << endl;
                     Pout << "nuf = " << nuf << endl;
                     Pout << "voidfraction = " << voidfraction << endl;
                     Pout << "Rep = " << Rep << endl;
-                    Pout << "beta = " << beta << endl;
+                    Pout << "betaP = " << betaP << endl;
                     Pout << "drag = " << drag << endl;
                 }
 
@@ -221,7 +249,7 @@ void GidaspowDrag::setForce() const
                     vValues.append(drag);   //first entry must the be the force
                     vValues.append(Ur);
                     sValues.append(Rep);
-                    sValues.append(beta);
+                    sValues.append(betaP);
                     sValues.append(voidfraction);
                     particleCloud_.probeM().writeProbe(index, sValues, vValues);
                 }
@@ -230,7 +258,21 @@ void GidaspowDrag::setForce() const
             // set force on particle
             if(treatExplicit_) for(int j=0;j<3;j++) expForces()[index][j] += drag[j];
             else  for(int j=0;j<3;j++) impForces()[index][j] += drag[j];
-            for(int j=0;j<3;j++) DEMForces()[index][j] += drag[j];
+
+            // set Cd
+            if(implDEM_)
+            {
+                for(int j=0;j<3;j++) fluidVel()[index][j]=Ufluid[j];
+
+                if (modelType_=="B" && cellI > -1)
+                    Cds()[index][0] = Vs*betaP/voidfraction*scaleDrag_;
+                else
+                    Cds()[index][0] = Vs*betaP*scaleDrag_;
+
+            }else{
+                for(int j=0;j<3;j++) DEMForces()[index][j] += drag[j];
+            }
+
         //}// end if mask
     }// end loop particles
 }

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

@@ -80,6 +80,10 @@ private:
 
     bool interpolation_; // use interpolated field values
 
+    mutable scalar scaleDia_;
+
+    mutable scalar scaleDrag_;
+
 public:
 
     //- Runtime type information

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

@@ -115,8 +115,8 @@ void KochHillDrag::setForce() const
 {
     if (scaleDia_ > 1)
         Info << "KochHill using scale = " << scaleDia_ << endl;
-    else if (cg() > 1){
+    else if (particleCloud_.cg() > 1){
-        scaleDia_=cg();
+        scaleDia_=particleCloud_.cg();
         Info << "KochHill using scale from liggghts cg = " << scaleDia_ << endl;
     }
 
@@ -158,6 +158,7 @@ void KochHillDrag::setForce() const
             betaP = 0;
             Vs = 0;
             Ufluid =vector(0,0,0);
+            voidfraction=0;
 
             if (cellI > -1) // particle Found
             {
@@ -233,6 +234,7 @@ void KochHillDrag::setForce() const
                     Pout << "nuf = " << nuf << endl;
                     Pout << "voidfraction = " << voidfraction << endl;
                     Pout << "Rep = " << Rep << endl;
+                    Pout << "betaP = " << betaP << endl;
                     Pout << "drag = " << drag << endl;
                 }
 
@@ -257,10 +259,10 @@ void KochHillDrag::setForce() const
             {
                 for(int j=0;j<3;j++) fluidVel()[index][j]=Ufluid[j];
 
-                if (modelType_=="B")
+                if (modelType_=="B" && cellI > -1)
-                    Cds()[index][0] = Vs*betaP/voidfraction;
+                    Cds()[index][0] = Vs*betaP/voidfraction*scaleDrag_;
                 else
-                    Cds()[index][0] = Vs*betaP;
+                    Cds()[index][0] = Vs*betaP*scaleDrag_;
 
             }else{
                 for(int j=0;j<3;j++) DEMForces()[index][j] += drag[j];

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

@@ -140,12 +140,10 @@ void KochHillRWDrag::setForce() const
     // realloc the arrays
     reAllocArrays();
 
-Info << "huhu" << endl;
-
     if (scale_ > 1)
         Info << "KochHillRW using scale = " << scale_ << endl;
-    else if (cg() > 1){
+    else if (particleCloud_.cg() > 1){
-        scale_=cg();
+        scale_=particleCloud_.cg();
         Info << "KochHillRW using scale from liggghts cg = " << scale_ << endl;
     }
 

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

@@ -99,7 +99,7 @@ void checkCouplingInterval::setForce() const
 
             if (cellI > -1) // particle Found
             {
-                scaledRad = particleCloud_.radius(index)/cg();
+                scaledRad = particleCloud_.radius(index)/particleCloud_.cg();
                 tauP = rhoP_*4*scaledRad*scaledRad/
                         (18 * nufField[cellI] * rho_[cellI]);
                 minTauP = min(minTauP,tauP);

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

@@ -74,9 +74,8 @@ fieldTimeAverage::fieldTimeAverage
     // create time average scalar fields
     scalarFields_.setSize(scalarFieldNames_.size());
 
-    if (propsDict_.found("startTime")){
+    if (propsDict_.found("startTime"))
         startTime_=readScalar(propsDict_.lookup("startTime"));
-    }
 
     for (int i=0;i < scalarFieldNames_.size(); i++)
     {

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

@@ -73,8 +73,7 @@ private:
     mutable PtrList<volVectorField> vectorFields_;
 
     mutable double nrAverages_;
-
+    
-
 public:
 
     //- Runtime type information

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

@@ -85,8 +85,8 @@ forceModel::forceModel
     ),
     coupleForce_(true),
     modelType_(sm.modelType()),
-    cg_(1.),
+    probeIt_(sm.probeM().active()),
-    probeIt_(sm.probeM().active())
+    requiresEx_(false)
 {}
 
 
@@ -130,6 +130,21 @@ void forceModel::manipulateScalarField(volScalarField& field) const
     Info << "no scalar manipulation done" << endl;
     // do nothing
 }
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void forceModel::repartitionImExForces() const
+{
+  if(particleCloud_.imExSplitFactor()<1.0)
+  {
+    Info << "Will re-partition split of implicit and explicit forces: alpha = " 
+         << particleCloud_.imExSplitFactor() << endl;
+    // Update implicit particle
+    expParticleForces_ += (1.0-particleCloud_.imExSplitFactor())*impParticleForces_;
+    impParticleForces_ *= particleCloud_.imExSplitFactor();
+  }
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -75,9 +75,9 @@ protected:
 
         const word modelType_;
 
-        mutable scalar cg_;
-
         bool probeIt_;
+        
+        bool requiresEx_;           //requires a orientation vector
 
 public:
 
@@ -150,17 +150,12 @@ public:
 
         inline const bool& coupleForce() const { return coupleForce_;};
 
-        inline const scalar cg() const { return cg_; };
-
-        void const setCG(double cg) const
-        { 
-            cg_ = cg;
-            Info << "cg is set to: " << cg_ << endl;
-        };
-
         void const setImpDEMdrag() const {particleCloud_.impDEMdrag_=true;};
-
+        
-
+        virtual inline bool& requiresEx() { return requiresEx_;};
+        
+    //Repartition Implixit/Explicit forces
+        void repartitionImExForces() const;
 };
 
 

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

@@ -100,7 +100,15 @@ void noDrag::setForce() const
                 if(treatExplicit_) for(int j=0;j<3;j++) expForces()[index][j] = 0.;
                 else  for(int j=0;j<3;j++) impForces()[index][j] = 0.;
             }
-            if(noDEMForce_)for(int j=0;j<3;j++) DEMForces()[index][j] = 0.;
+            if(noDEMForce_)
+            {
+                for(int j=0;j<3;j++) DEMForces()[index][j] = 0.;
+                if(particleCloud_.impDEMdrag())
+                {
+                    Cds()[index][0] = 0.;
+                    for(int j=0;j<3;j++) fluidVel()[index][j] = 0.;
+                }
+            }
         }        
     }
 }

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.C

@@ -89,73 +89,9 @@ execute::execute
     // check if verbose
     if (propsDict_.found("verbose")) verbose_=true;
 
-    bool addBlank = true;  // std no blanks after each word
+    parseCommandList(commandList_, labelList_, scalarList_, command_, propsDict_, timeStamp_);
-    fileName add;
-    label numberCount=0;   // nr of scalars inserted to command
-    label labelCount=0;   // nr of labels inserted to command
-
-    forAll(commandList_,i)
-    {
-        add = word(commandList_[i]);
-
-        //- handle symbols
-        if (add == "$couplingInterval")
-        {
-            char h[50];
-            sprintf(h,"%d",particleCloud_.dataExchangeM().couplingInterval());
-            add = h;
-        }
-        else if (add=="dot")    add = ".";
-        else if (add=="dotdot") add = "..";
-        else if (add=="slash")  add = "/";
-        else if (add=="noBlanks")  // no blanks after the following words
-        {
-            add = "";
-            addBlank = false;
-        }else if (add=="blanks") // add a blank here and after the following words
-        {
-            add = "";
-            addBlank = true;
-        }else if (add=="timeStamp") // next command will be a number read from scalarList_
-        {
-            add = "";
-            timeStamp_=true;
-        }else if (add=="number") // next command will be a number read from scalarList_
-        {
-            if (!propsDict_.found("scalars"))
-            {
-                FatalError<<"you want to use a number in the command\n - specify a scalar list with all numbers"
-                << abort(FatalError);
-
-            }
-            char h[50];
-            sprintf(h,"%f",scalarList_[numberCount]);
-            add = h;
-            numberCount ++;
-        }else if (add=="label") // next command will be a number read from scalarList_
-        {
-            if (!propsDict_.found("labels"))
-            {
-                FatalError<<"you want to use a label in the command\n - specify a label list with all numbers"
-                << abort(FatalError);
-
-            }
-            char h[50];
-            sprintf(h,"%d",labelList_[labelCount]);
-            add = h;
-            labelCount ++;
-        }
-
-        // compose command
-        if (addBlank)
-        {
-            command_ += add + " ";
-        }else
-        {
-            command_ += add;
-        }
-    }
     Info << "liggghtsCommand " << command_ << endl;
+
     strCommand_=string(command_);
 
     checkTimeMode(propsDict_);
@@ -172,7 +108,7 @@ execute::~execute()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const char* execute::command()
+const char* execute::command(int commandLine)
 {
     return strCommand_.c_str();
 }

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/execute/execute.H

@@ -98,7 +98,7 @@ public:
 
     word name(){return myName_;};
 
-    const char* command();
+    const char* command(int);
 
     bool runCommand(int);
 

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C

@@ -71,6 +71,7 @@ liggghtsCommandModel::liggghtsCommandModel
     lastCouplingStep_(-1),
     couplingStepInterval_(0),
     exactTiming_(false),
+    commandLines_(1),
     verbose_(false)
 {}
 
@@ -236,6 +237,74 @@ DynamicList<scalar> liggghtsCommandModel::executionsWithinPeriod(scalar TSstart,
 
     return executions;
 }
+
+void liggghtsCommandModel::parseCommandList(wordList& commandList,labelList& labelList,scalarList& scalarList,word& command, dictionary& propsDict, bool timeStamp)
+{
+    bool addBlank = true;  // std no blanks after each word
+    fileName add;
+    label numberCount=0;   // nr of scalars inserted to command
+    label labelCount=0;   // nr of labels inserted to command
+
+    forAll(commandList,i)
+    {
+        add = word(commandList[i]);
+
+        //- handle symbols
+        if (add == "$couplingInterval")
+        {
+            char h[50];
+            sprintf(h,"%d",particleCloud_.dataExchangeM().couplingInterval());
+            add = h;
+        }
+        else if (add=="dot")    add = ".";
+        else if (add=="dotdot") add = "..";
+        else if (add=="slash")  add = "/";
+        else if (add=="noBlanks")  // no blanks after the following words
+        {
+            add = "";
+            addBlank = false;
+        }else if (add=="blanks") // add a blank here and after the following words
+        {
+            add = "";
+            addBlank = true;
+        }else if (add=="timeStamp") // next command will be a number read from labelList
+        {
+            add = "";
+            timeStamp=true;
+        }else if (add=="number") // next command will be a number read from labelList
+        {
+            /*if (!propsDict.found("scalars"))
+            {
+                FatalError<<"you want to use a number in the command\n - specify a scalar list with all numbers"
+                << abort(FatalError);
+            }*/
+            char h[50];
+            sprintf(h,"%f",scalarList[numberCount]);
+            add = h;
+            numberCount ++;
+        }else if (add=="label") // next command will be a number read from labelList
+        {
+            /*if (!propsDict.found("labels"))
+            {
+                FatalError<<"you want to use a label in the command\n - specify a label list with all numbers"
+                << abort(FatalError);
+            }*/
+            char h[50];
+            sprintf(h,"%d",labelList[labelCount]);
+            add = h;
+            labelCount ++;
+        }
+
+        // compose command
+        if (addBlank)
+        {
+            command += add + " ";
+        }else
+        {
+            command += add;
+        }
+    }
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.H

@@ -89,6 +89,8 @@ protected:
 
     bool exactTiming_;
 
+    label commandLines_;
+
     bool verbose_;
 
 public:
@@ -142,7 +144,7 @@ public:
 
     // Member Functions
 
-        virtual const char* command()=0;
+        virtual const char* command(int)=0;
 
         void checkTimeMode(dictionary&);
 
@@ -168,6 +170,9 @@ public:
 
         bool exactTiming(){return exactTiming_;};
 
+        label commandLines(){return commandLines_;};
+
+        void parseCommandList(wordList&, labelList&, scalarList&, word&, dictionary&, bool=false);
 };
 
 

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.C

@@ -112,7 +112,7 @@ readLiggghtsData::~readLiggghtsData()
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-const char* readLiggghtsData::command()
+const char* readLiggghtsData::command(int commandLine)
 {
     char h[50];
     sprintf(h,"_%d",insertionNr_);

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/readLiggghtsData/readLiggghtsData.H

@@ -93,7 +93,7 @@ public:
 
     word name(){return myName_;};
 
-    const char* command();
+    const char* command(int);
 
     bool runCommand(int);
 

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.C

@@ -93,7 +93,7 @@ runLiggghts::~runLiggghts()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const char* runLiggghts::command()
+const char* runLiggghts::command(int commandLine)
 {
     return strCommand_.c_str();
 }

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/runLiggghts/runLiggghts.H

@@ -88,7 +88,7 @@ public:
 
     // Member Functions
 
-    const char* command();
+    const char* command(int);
 
     string createCommand(word="",int=0,word="",word="",word="",word="");
 

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.C

@@ -115,7 +115,7 @@ writeLiggghts::~writeLiggghts()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const char* writeLiggghts::command()
+const char* writeLiggghts::command(int commandLine)
 {
     return strCommand_.c_str();
 }

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.H

@@ -93,7 +93,7 @@ public:
 
     // Member Functions
 
-    const char* command();
+    const char* command(int);
 
     bool runCommand(int);
 

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U

@@ -51,29 +51,12 @@ boundaryField
         value           uniform (0 0 0);*/
 
         // 2.1.x
+		//type              uniformFixedValue;
         type            uniformFixedValueVoidfraction;
         uniformValue     table
         (
             (0.000 (0 0 0.002))
-            (0.020 (0 0 0.002))
+            (0.100 (0 0 0.030))
-            (0.021 (0 0 0.004))
-            (0.040 (0 0 0.004))
-            (0.041 (0 0 0.006))
-            (0.060 (0 0 0.006))
-            (0.061 (0 0 0.008))
-            (0.080 (0 0 0.008))
-            (0.081 (0 0 0.010))
-            (0.100 (0 0 0.010))
-            (0.101 (0 0 0.012))
-            (0.120 (0 0 0.012))
-            (0.121 (0 0 0.014))
-            (0.140 (0 0 0.014))
-            (0.141 (0 0 0.016))
-            (0.160 (0 0 0.016))
-            (0.161 (0 0 0.018))
-            (0.180 (0 0 0.018))
-            (0.181 (0 0 0.020))
-            (0.200 (0 0 0.020))
         );
     }
 

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/p

@@ -35,8 +35,9 @@ boundaryField
     outlet
     {
         //type            zeroGradient;
+
         type            fixedValue;
-        value           uniform 100000;
+        value           uniform 1.0e5;
     }
 }
 

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/rho

@@ -16,7 +16,7 @@ FoamFile
 
 dimensions      [1 -3 0 0 0 0 0];
 
-internalField   uniform 30;
+internalField   uniform 10;
 
 boundaryField
 {

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/voidfraction

@@ -28,10 +28,14 @@ boundaryField
     outlet
     {
         type            zeroGradient;
+        //type            fixedValue;
+        //value           uniform 1;
     }
     inlet
     {
         type            zeroGradient;
+        //type            fixedValue;
+        //value           uniform 1;
     }
 }
 

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties

@@ -25,42 +25,44 @@ FoamFile
 
 //===========================================================================//
 // sub-models & settings
+solveFlow true;
 
 modelType "A"; // A or B
 
 couplingInterval 100;
 
-voidFractionModel divided;
+voidFractionModel divided;//centre;//
 
-locateModel engine;
+locateModel engine;//turboEngineM2M;//
 
 meshMotionModel noMeshMotion;
 
 regionModel allRegion;
 
-IOModel "basicIO";
+IOModel basicIO;
 
 probeModel off;
 
-dataExchangeModel twoWayMPI;//twoWayFiles;//oneWayVTK;//
+dataExchangeModel twoWayMPI;//twoWayM2M;//twoWayFiles;//oneWayVTK;//
 
 averagingModel dense;//dilute;//
 
 clockModel off;//standardClock;//
 
-smoothingModel off;
+smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; //
 
 forceModels
 (
-    //GidaspowDrag
+    //KochHillDrag
-    //DiFeliceDrag
-    KochHillDrag
     //BeetstraDrag
+    //noDrag
     gradPForce
     viscForce
+    GidaspowDrag
     //Archimedes
     //volWeightedAverage
     //totalMomentumExchange
+    particleCellVolume
 );
 
 momCoupleModels
@@ -73,6 +75,21 @@ turbulenceModelType "RASProperties";//"LESProperties";//
 //===========================================================================//
 // sub-model properties
 
+localPSizeDiffSmoothingProps
+{
+    lowerLimit 0.1;
+    upperLimit 1e10;
+    dSmoothingLength 1.5e-3;
+    Csmoothing    1.0; 
+}
+
+constDiffSmoothingProps
+{
+    lowerLimit 0.1;
+    upperLimit 1e10;
+    smoothingLength 1.5e-3;
+}
+
 implicitCoupleProps
 {
     velFieldName "U";
@@ -91,7 +108,7 @@ gradPForceProps
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
     velocityFieldName "U";
-    //interpolation;
+    interpolation;
 }
 
 viscForceProps
@@ -100,6 +117,11 @@ viscForceProps
     densityFieldName "rho";
     interpolation;
 }
+noDragProps
+{
+    //noDEMForce;
+    //keepCFDForce;
+}
 volWeightedAverageProps
 {
     scalarFieldNames
@@ -123,22 +145,24 @@ totalMomentumExchangeProps
 }
 GidaspowDragProps
 {
-    velFieldName "U";
+    verbose;
-    densityFieldName "rho";
-}
-DiFeliceDragProps
-{
     velFieldName "U";
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    phi 1;
+    //interpolation; // this case does not like interpolation
+    implDEM_;
 }
 
 KochHillDragProps
 {
+    //verbose;
     velFieldName "U";
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
-    //verbose;
+    //interpolation;
+    implDEM_;
+    verbose;
 }
 
 BeetstraDragProps
@@ -162,6 +186,13 @@ virtualMassForceProps
     densityFieldName "rho";
 }
 
+particleCellVolumeProps
+{
+    upperThreshold 0.999;
+    lowerThreshold 0.;
+    verbose;
+}
+
 oneWayVTKProps
 {
     couplingFilename "vtk_out%4.4d.vtk";
@@ -175,7 +206,7 @@ twoWayFilesProps
 
 centreProps
 {
-    alphaMin 0.10;
+    alphaMin 0.1;
 }
 
 engineProps
@@ -185,7 +216,7 @@ engineProps
 
 dividedProps
 {
-    alphaMin 0.1;
+    alphaMin 0.05;
     scaleUpVol 1.0;
 }
 

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

@@ -5,7 +5,7 @@ clc;
 %====================================%
 % simulation data 1
 %====================================%
-rhoG = 30			% density in kg/m3
+rhoG = 10			% density in kg/m3
 %path = '../probes/0/p';
 path = '../postProcessing/probes/0/p';
 columns=22;
@@ -13,7 +13,7 @@ headerlines=4;
 data = loaddata(path,columns,headerlines);
 data=transpose(data);
 [x,y]=size(data)
-dp_sim = (data(:,2)-data(:,y))*rhoG;
+dp_sim = (data(:,2)-data(:,y))*rhoG; %conversion to Pa!
 t_sim = data(:,1);
 %fprintf('final pressureDrop of sim = %f Pa\n',dp_sim(length(dp_sim)) )
 
@@ -25,18 +25,18 @@ t_sim = data(:,1);
 % Ergun Equation
 %===================
 fprintf('\ncalc Ergun eqn:\n')
-dp = 0.0008			% particle diameter
+dp = 0.0005			% particle diameter
 phip = 1			% sphericity
-epsilon = 0.451335              % void fraction
+epsilon = 1-0.518497 % void fraction NOTE: varies with voidfraction!!!
 Ustart = 0.002
-Uend = 0.02
+Uend = 0.03
 timeStepSize = 0.001;            % time interval of pressure data
 Tstart = 0;
 Tend = t_sim(length(t_sim));
 deltaU=(Uend-Ustart)/((Tend-Tstart)/timeStepSize);
 U = Ustart+deltaU:deltaU:Uend;  % velocity over time
 Ua = U / epsilon;		% physical velocity
-L = 0.0156			% length of bed
+L = 0.0135			% length of bed
 nuG = 1.5*10^-4			% kinemat Visk in m2/s
 muG = nuG*rhoG			% dynam visc in Pa s
 
@@ -55,14 +55,31 @@ fprintf('so the result does not depend on density\n')
 %==================================
 rhoP = 7000                      % particle density in kg/m3
 g = 9.81                        % gravity m/s2
+
+% WARNING!!!
+% these equations can give wrong results if applied to the wrong regime!!!
+
 Umf = dp^2*(rhoP-rhoG)*g/(150*muG)*(epsilon^3*phip^2)/(1-epsilon)
-ReMF = Umf*dp*rhoG/muG 			% must be <20 !!!
+ReMF = Umf*dp*rhoG/muG
-%Umf = sqrt(phip*dp^2/1.75*(rhoP-rhoG)/rhoG*g*epsilon^3) % Re>1000
+#Umf = sqrt(dp*(rhoP-rhoG)*g/(1.75*rhoG)*epsilon^3*phip)
+#ReMF = Umf*dp*rhoG/muG
+
+
+if(ReMF<20)
+    fprintf('applying eqn1 for Umf.\n')
+elseif(ReMF>20 && ReMF<1000)
+    fprintf('applying eqn1 for Umf.\n')
+elseif (ReMF>=1000)
+    fprintf('applying eqn2 for Umf.\n')
+    Umf = sqrt(dp*(rhoP-rhoG)*g/(1.75*rhoG)*epsilon^3*phip);
+    ReMF = Umf*dp*rhoG/muG;
+end
 
 dpUmf= L * (
                 150*((1-epsilon)^2/epsilon^3)*((muG.*Umf)/(phip*dp)^2) 
               +1.75*((1-epsilon)/epsilon^3)*((rhoG.*Umf.^2)/(phip*dp))
-        );
+        )
+dpUmf=(L*(1-epsilon)*(rhoP-rhoG)*g)
 
 %====================================%
 % plot data
@@ -80,7 +97,7 @@ figure(1)
 plot(U,dpErgun,U,dp_sim,[Umf,Uend],dpUmf*ones(1,2))
 title("Ergun pressure drop vs. simulation")
 a=strcat("analytical (Ergun), Umf=",num2str(Umf),", dpUmf=",num2str(dpUmf));
-legend(a,"simulation")
+legend(a,"simulation","analyt. deltaP at Umf")
 xlabel("velocity in [m/s]")
 ylabel("pressure drop [Pa]")
 axis([0,Uend,0,dpErgun(length(dpErgun))])

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict

@@ -23,13 +23,13 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.2;//0.01;
+endTime         0.1;
 
 deltaT          0.001;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.005;//0.01;
+writeInterval   0.01;
 
 purgeWrite      0;
 
@@ -54,6 +54,7 @@ libs ( "libfiniteVolumeCFDEM.so" );
 
 functions
 (
+
     probes
     {
         type        probes;
@@ -87,7 +88,7 @@ functions
         );
 
         // Fields to be probed
-        fields ( p U voidfraction volAverage_voidfraction);
+        fields ( p U voidfraction volAverage_voidfraction voidfractionNext voidfractionPrev);
 
         // Write at same frequency as fields
         outputControl   timeStep;//outputTime;

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/decomposeParDict

@@ -30,7 +30,7 @@ numberOfSubdomains  4;
 
 simpleCoeffs
 {
-    n           (1 2 2);
+    n           (2 2 1);
     delta       0.001;
 }
 

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    p
+    "(p)"
     {
         solver          PCG;
         preconditioner  DIC;
@@ -33,7 +33,7 @@ solvers
         relTol          0;
     }
 
-    U
+    "(U|k|epsilon|R|nuTilda)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -41,34 +41,10 @@ solvers
         relTol          0;
     }
 
-    k
+    "(voidfraction|Us|Ksl|dSmoothing)"
     {
-        solver          PBiCG;
+        solver          PCG;
-        preconditioner  DILU;
+        preconditioner  DIC;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    epsilon
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    R
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    nuTilda
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
         tolerance       1e-05;
         relTol          0;
     }

tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_init

@@ -7,9 +7,9 @@ echo both
 #######################################################
 # variables                                           #
 #######################################################
-variable cg equal 2
+variable cg equal 1
-variable r0 equal 0.0004
+variable r0 equal 0.00025
-variable nPorg equal 19531
+variable nPorg equal 71000
 variable nPscal equal ${nPorg}/(${cg}*${cg}*${cg})
 #######################################################
 
@@ -22,6 +22,7 @@ boundary	m m m
 newton		off
 
 units		si
+processors	2 2 1
 
 region		reg block -0.015 0.015 -0.015 0.015 -0.001 0.0554 units box
 create_box	1 reg
@@ -36,24 +37,23 @@ fix 		m1 all property/global youngsModulus peratomtype 5.e6
 fix 		m2 all property/global poissonsRatio peratomtype 0.45
 fix 		m3 all property/global coefficientRestitution peratomtypepair 1 0.3
 fix 		m4 all property/global coefficientFriction peratomtypepair 1 0.5
-fix 		m5 all property/global characteristicVelocity scalar 1.
 
 #pair style
-pair_style gran model hooke tangential history
+pair_style gran model hertz tangential history #Hertzian without cohesion
 pair_coeff	* *
 
 #timestep, gravity
-timestep	0.000005
+timestep	0.00001
 fix		gravi all gravity 9.81 vector 0.0 0.0 -1.0
 
 #walls
-fix zwalls1 all wall/gran model hooke tangential history primitive type 1  zplane 0.0
+fix zwalls1 all wall/gran model hertz tangential history primitive type 1  zplane 0.0
-fix zwalls2 all wall/gran model hooke tangential history primitive type 1  zplane 0.0553
+fix zwalls2 all wall/gran model hertz tangential history primitive type 1  zplane 0.0553
-fix cylwalls all wall/gran model hooke tangential history primitive type 1 zcylinder 0.01385 0. 0.
+fix cylwalls all wall/gran model hertz tangential history primitive type 1 zcylinder 0.01385 0. 0.
 
 #particle distributions and insertion
 region		bc cylinder z 0.0 0.0 0.012 0. 0.055 units box
-fix		    pts1 all particletemplate/sphere 1 atom_type 1 density constant 200 radius constant ${r0}
+fix		    pts1 all particletemplate/sphere 1 atom_type 1 density constant 7000 radius constant ${r0}
 fix		    pdd1 all particledistribution/discrete 1.  1 pts1 1.0
 
 fix		ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -1. insert_every once overlapcheck yes all_in yes particles_in_region ${nPscal} region bc

tutorials/cfdemSolverPiso/ErgunTestCG/DEM/in.liggghts_resume

@@ -7,9 +7,9 @@ echo both
 #######################################################
 # variables                                           #
 #######################################################
-variable cg equal 2
+variable cg equal 3
-variable r0 equal 0.0004
+variable r0 equal 0.00025
-variable nPorg equal 19531
+variable nPorg equal 71000
 variable nPscal equal ${nPorg}/(${cg}*${cg}*${cg})
 #######################################################
 
@@ -22,28 +22,24 @@ boundary	m m m
 newton		off
 
 units		si
-processors	1 2 2
+processors	2 2 1
 
 #read the restart file
 read_restart 	../DEM/liggghts.restart_coarseGraining_${cg}
 
-#do not do this here, the simulation box is in the restart file!
+neighbor	0.001 bin
-#region		reg block -0.015 0.015 -0.015 0.015 -0.001 0.0554 units box
-#create_box	1 reg
-
-neighbor	0.0005 bin
 neigh_modify	delay 0
 
+
 #Material properties required for new pair styles
 
 fix 		m1 all property/global youngsModulus peratomtype 5.e6
 fix 		m2 all property/global poissonsRatio peratomtype 0.45
 fix 		m3 all property/global coefficientRestitution peratomtypepair 1 0.3
 fix 		m4 all property/global coefficientFriction peratomtypepair 1 0.5
-fix 		m5 all property/global characteristicVelocity scalar 1.
 
 #pair style
-pair_style gran model hooke tangential history #Hertzian without cohesion
+pair_style gran model hertz tangential history #Hertzian without cohesion
 pair_coeff	* *
 
 #timestep, gravity
@@ -51,16 +47,14 @@ timestep	0.00001
 fix		gravi all gravity 9.81 vector 0.0 0.0 -1.0
 
 #walls
-fix zwalls1 all wall/gran model hooke tangential history primitive type 1  zplane 0.0
+fix zwalls1 all wall/gran model hertz tangential history primitive type 1  zplane 0.0
-fix zwalls2 all wall/gran model hooke tangential history primitive type 1  zplane 0.0553
+fix zwalls2 all wall/gran model hertz tangential history primitive type 1  zplane 0.0553
-fix cylwalls all wall/gran model hooke tangential history primitive type 1 zcylinder 0.01385 0. 0.
+fix cylwalls all wall/gran model hertz tangential history primitive type 1 zcylinder 0.01385 0. 0.
 
-# change the particles density
-set             group all density 7000
 
 #cfd coupling
 fix		cfd all couple/cfd couple_every 100 mpi
-fix		cfd2 all couple/cfd/force
+fix		cfd2 all couple/cfd/force #/implicit
 
 #apply nve integration to all particles that are inserted as single particles
 fix		integr all nve/sphere
@@ -69,11 +63,11 @@ fix		integr all nve/sphere
 compute         centerOfMass all com
 
 #compute total dragforce
-compute		dragtotal all reduce sum f_dragforce[1] f_dragforce[2] f_dragforce[3]
+#compute		dragtotal all reduce sum f_dragforce[1] f_dragforce[2] f_dragforce[3]
 
 #screen output
 compute		1 all erotate/sphere
-thermo_style	custom step atoms ke c_1 vol c_centerOfMass[3] c_dragtotal[1] c_dragtotal[2] c_dragtotal[3]
+thermo_style	custom step atoms ke c_1 vol c_centerOfMass[3] #c_dragtotal[1] c_dragtotal[2] c_dragtotal[3]
 thermo		10
 thermo_modify	lost ignore norm no
 compute_modify	thermo_temp dynamic yes
@@ -81,7 +75,7 @@ compute_modify	thermo_temp dynamic yes
 #insert the first particles so that dump is not empty
 dump 		myDump all stl 1 post/dump_*.stl
 #run		1
-dump		dmp all custom 5000 ../DEM/post/dump*.liggghts_restart id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz f_dragforce[1] f_dragforce[2] f_dragforce[3] radius 
+dump		dmp all custom 1000 ../DEM/post/dump*.liggghts_restart id type type x y z vx vy vz fx fy fz radius #f_dragforce[1] f_dragforce[2] f_dragforce[3] 
 undump		myDump
 
 run		1

tutorials/cfdemSolverPiso/ErgunTestCG/DEM/post/dummy

@@ -1 +0,0 @@
-dummyfile

tutorials/cfdemSolverPiso/ErgunTestCG/parCFDDEMrun.sh

@@ -21,7 +21,7 @@ logfileName="log_$headerText"
 solverName="cfdemSolverPiso"
 nrProcs="4"
 machineFileName="none"   # yourMachinefileName | none
-debugMode="off"          # on | off | strict
+debugMode="off"          # on | off| strict
 testHarnessPath="$CFDEM_TEST_HARNESS_PATH"
 runOctave="true"
 postproc="false"
@@ -82,21 +82,10 @@ fi
 
 #- clean up case
 echo "deleting data at: $casePath :\n"
-rm -r $casePath/CFD/0.*
+source $WM_PROJECT_DIR/bin/tools/CleanFunctions
-rm -r $casePath/CFD/log.*
+cd $casePath/CFD
-rm -r $casePath/CFD/octave/octave-core
+cleanCase
-rm -r $casePath/CFD/VTK
-rm -r $casePath/CFD/processor*
-rm -r $casePath/CFD/couplingFiles/*
-rm -r $casePath/DEM/post/*
-rm -r $casePath/DEM/log.*
-rm -r $casePath/DEM/liggghts.restartCFDEM*
-rm -r $casePath/CFD/probes*
-rm -r $casePath/CFD/postProcessing
-rm -r $casePath/CFD/particles
-rm -r $casePath/CFD/lagrangian
 rm -r $casePath/CFD/clockData
+rm -r $casePath/DEM/post/*
+(cd $casePath/DEM/post && touch dummy)
 echo "done"
-
-#- preserve post directory
-echo "dummyfile" >> $casePath/DEM/post/dummy

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/U

@@ -55,24 +55,6 @@ boundaryField
         uniformValue     table
         (
             (0.000 (0 0 0.002))
-            /*(0.010 (0 0 0.002))
-            (0.011 (0 0 0.004))
-            (0.020 (0 0 0.004))
-            (0.021 (0 0 0.006))
-            (0.030 (0 0 0.006))
-            (0.031 (0 0 0.008))
-            (0.040 (0 0 0.008))
-            (0.041 (0 0 0.010))
-            (0.050 (0 0 0.010))
-            (0.051 (0 0 0.012))
-            (0.060 (0 0 0.012))
-            (0.061 (0 0 0.014))
-            (0.070 (0 0 0.014))
-            (0.071 (0 0 0.016))
-            (0.080 (0 0 0.016))
-            (0.081 (0 0 0.018))
-            (0.090 (0 0 0.018))
-            (0.091 (0 0 0.020))*/
             (0.100 (0 0 0.020))
         );
     }

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/couplingProperties

@@ -62,6 +62,7 @@ forceModels
     //volWeightedAverage
     //totalMomentumExchange
     //particleCellVolume
+    //fieldTimeAverage
 );
 
 momCoupleModels
@@ -107,7 +108,7 @@ gradPForceProps
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
     velocityFieldName "U";
-    //interpolation;
+    interpolation;
 }
 
 viscForceProps
@@ -143,6 +144,7 @@ GidaspowDragProps
     velFieldName "U";
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    interpolation;
     phi 1;
 }
 DiFeliceDragProps
@@ -158,6 +160,7 @@ KochHillDragProps
     velFieldName "U";
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    interpolation;
 }
 
 BeetstraDragProps
@@ -201,7 +204,7 @@ twoWayFilesProps
 
 centreProps
 {
-    alphaMin 0.10;
+    alphaMin 0.1;
 }
 
 engineProps

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/liggghtsCommands

@@ -25,7 +25,7 @@ FoamFile
 
 liggghtsCommandModels
 (
-   runLiggghts
+    runLiggghts
 );
 // ************************************************************************* //
 

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

@@ -30,7 +30,7 @@ phip = 1			% sphericity
 epsilon = 0.451335              % void fraction
 Ustart = 0.002
 Uend = 0.02
-timeStepSize = 0.001;            % time interval of pressure data
+timeStepSize = 0.0005;            % time interval of pressure data
 Tstart = 0;
 Tend = t_sim(length(t_sim));
 deltaU=(Uend-Ustart)/((Tend-Tstart)/timeStepSize);

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/controlDict

@@ -25,7 +25,7 @@ stopAt          endTime;
 
 endTime         0.1;
 
-deltaT          0.001;
+deltaT          0.0005;
 
 writeControl    adjustableRunTime;
 

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/fvSolution

@@ -41,7 +41,7 @@ solvers
         relTol          0;
     }
 
-    "(voidfraction|Us|Ksl|dSmoothing)"
+    "(voidfraction|Us|Ksl|dSmoothing|UsNext|voidfractionNext)"
     {
         solver          PCG;
         preconditioner  DIC;

tutorials/cfdemSolverPiso/ErgunTestMPI/DEM/in.liggghts_resume

@@ -65,7 +65,7 @@ compute_modify	thermo_temp dynamic yes
 #insert the first particles so that dump is not empty
 dump 		myDump all stl 1 post/dump_*.stl
 #run		1
-dump		dmp all custom 5000 ../DEM/post/dump*.liggghts_restart id type type x y z vx vy vz fx fy fz f_dragforce[1] f_dragforce[2] f_dragforce[3] radius 
+dump		dmp all custom 5000 ../DEM/post/dump*.liggghts_restart id type type x y z vx vy vz fx fy fz f_dragforce[1] f_dragforce[2] f_dragforce[3] radius
 undump		myDump
 
 run		1

tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/couplingProperties

@@ -130,6 +130,7 @@ DiFeliceDragProps
     velFieldName "U";
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
 }
 
 KochHillDragProps

tutorials/cfdemSolverPiso/settlingTestMPI/CFD/constant/couplingProperties

@@ -93,6 +93,8 @@ DiFeliceDragProps
     velFieldName "U";
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
+    verbose;
 }
 
 SchillerNaumannDragProps

tutorials/cfdemSolverPiso/settlingTestMPI/CFD/system/fvSolution

@@ -33,7 +33,7 @@ solvers
         relTol          0;
     }
 
-    U
+    "(U|k|epsilon|R|nuTilda)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -41,55 +41,13 @@ solvers
         relTol          0;
     }
 
-    k
+    "(voidfraction|Ksl|UsNext|voidfractionNext)"
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    epsilon
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    R
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    nuTilda
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    voidfraction
     {
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-09;
         relTol          1e-06;
     }
-
-    Ksl
-    {
-        $voidfraction
-    }
-
-    Us
-    {
-        $voidfraction
-    }
 }
 
 PISO