Merge pull request #80 from ParticulateFlow/release

Release 19.02

.gitignore

@@ -7,5 +7,6 @@ log.*
 *~
 
 **/linux*Gcc*/
+**/.vscode
 
 lnInclude

applications/solvers/cfdemSolverMultiphase/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wclean libso multiphaseMixture
+wclean
+
+#------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphase/Allwmake

@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Parse arguments for library compilation
+targetType=libso
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+set -x
+
+wmake $targetType multiphaseMixture
+wmake
+
+#------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphase/Make/files

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

applications/solvers/cfdemSolverMultiphase/Make/options

@@ -0,0 +1,30 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -ImultiphaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lcfdemMultiphaseInterFoam \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverMultiphase/UEqn.H

@@ -0,0 +1,61 @@
+const surfaceScalarField& rhoPhi(mixture.rhoPhi());
+
+volScalarField muEff = rho*(turbulence->nu() + turbulence->nut());
+
+if (modelType == "A")
+    muEff *= voidfraction;
+
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rhoEps, U) - fvm::Sp(fvc::ddt(rhoEps),U)
+  + fvm::div(rhoPhi, U) - fvm::Sp(fvc::div(rhoPhi),U)
+  //+ particleCloud.divVoidfractionTau(U, voidfraction)
+  - fvm::laplacian(muEff, U) - fvc::div(muEff*dev2(fvc::grad(U)().T()))
+  ==
+    fvOptions(rho, U)
+  - fvm::Sp(Ksl,U)
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve
+    (
+        UEqn
+     ==
+        fvc::reconstruct
+        (
+            (- ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh)) * mesh.magSf()
+        )
+     +
+        fvc::reconstruct
+        (
+            mixture.surfaceTensionForce() * mesh.magSf()
+        ) * voidfraction
+     + Ksl*Us
+    );
+
+    fvOptions.correct(U);
+}
+else if (pimple.momentumPredictor())
+{
+    solve
+    (
+        UEqn
+     ==
+        fvc::reconstruct
+        (
+            (
+                mixture.surfaceTensionForce()
+              - ghf*fvc::snGrad(rho)
+              - fvc::snGrad(p_rgh)
+            ) * mesh.magSf()
+        ) * voidfraction
+     + Ksl*Us
+    );
+
+    fvOptions.correct(U);
+}

applications/solvers/cfdemSolverMultiphase/additionalChecks.H

@@ -0,0 +1,17 @@
+// Additional solver-specific checks
+
+// Useful if one wants to e.g. initialize floating particles using the Archimedes model
+if (particleCloud.couplingProperties().found("unrestrictedForceModelSelection"))
+{
+    Warning << "Using unrestrictedForceModelSelection, results may be incorrect!" << endl;
+} else
+{
+    #include "checkModelType.H"
+}
+
+word modelType = particleCloud.modelType();
+
+if(!particleCloud.couplingProperties().found("useDDTvoidfraction"))
+{
+    Warning << "Suppressing ddt(voidfraction) is not recommended with this solver as it may generate incorrect results!" << endl;
+}

applications/solvers/cfdemSolverMultiphase/alphaCourantNo.H

@@ -0,0 +1,21 @@
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        mixture.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/cfdemSolverMultiphase.C

@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Application
+    cfdemSolverMultiphase
+
+Description
+    CFD-DEM solver for n incompressible fluids which captures the interfaces and
+    includes surface-tension and contact-angle effects for each phase. It is based
+    on the OpenFOAM(R)-4.x solver multiphaseInterFoam but extended to incorporate
+    DEM functionalities from the open-source DEM code LIGGGHTS.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseMixture.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+
+#include "cfdemCloud.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+
+    turbulence->validate();
+
+    // create cfdemCloud
+    cfdemCloud particleCloud(mesh);
+
+    #include "additionalChecks.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        if(particleCloud.solveFlow())
+        {
+            mixture.solve();
+            rho = mixture.rho();
+            rhoEps = rho * voidfraction;
+
+             // --- Pressure-velocity PIMPLE corrector loop
+            while (pimple.loop())
+            {
+                #include "UEqn.H"
+
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/correctPhi.H

@@ -0,0 +1,11 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"

applications/solvers/cfdemSolverMultiphase/createFields.H

@@ -0,0 +1,156 @@
+//===============================
+// particle interaction modelling
+//===============================
+
+Info<< "\nReading momentum exchange field Ksl\n" << endl;
+volScalarField Ksl
+(
+    IOobject
+    (
+        "Ksl",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+    //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
+);
+
+Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+volScalarField voidfraction
+(
+    IOobject
+    (
+        "voidfraction",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+);
+voidfraction.oldTime();
+
+Info<< "Reading particle velocity field Us\n" << endl;
+volVectorField Us
+(
+    IOobject
+    (
+        "Us",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+);
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+         mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(U*voidfraction) & mesh.Sf()
+);
+
+multiphaseMixture mixture(U, phi, voidfraction);
+
+// Need to store rho for ddt(rho, U)
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.rho()
+);
+rho.oldTime();
+
+volScalarField rhoEps ("rhoEps", rho * voidfraction);
+
+// Construct incompressible turbulence model
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, mixture)
+);
+
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    p_rgh + rho*gh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell
+(
+    p,
+    p_rgh,
+    pimple.dict(),
+    pRefCell,
+    pRefValue
+);
+
+if (p_rgh.needReference())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+}
+
+mesh.setFluxRequired(p_rgh.name());

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/Make/files

@@ -0,0 +1,5 @@
+phase/phase.C
+alphaContactAngle/alphaContactAngleFvPatchScalarField.C
+multiphaseMixture.C
+
+LIB = $(CFDEM_LIB_DIR)/libcfdemMultiphaseInterFoam

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/Make/options

@@ -0,0 +1,13 @@
+EXE_INC = \
+    -IalphaContactAngle \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+LIB_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
+(
+    Istream& is
+)
+:
+    theta0_(readScalar(is)),
+    uTheta_(readScalar(is)),
+    thetaA_(readScalar(is)),
+    thetaR_(readScalar(is))
+{}
+
+
+Istream& operator>>
+(
+    Istream& is,
+    alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
+    return is;
+}
+
+
+Ostream& operator<<
+(
+    Ostream& os,
+    const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    os  << tp.theta0_ << token::SPACE
+        << tp.uTheta_ << token::SPACE
+        << tp.thetaA_ << token::SPACE
+        << tp.thetaR_;
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(p, iF)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchScalarField(p, iF),
+    thetaProps_(dict.lookup("thetaProperties"))
+{
+    evaluate();
+}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, iF),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("thetaProperties")
+        << thetaProps_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphaContactAngleFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::alphaContactAngleFvPatchScalarField
+
+Description
+    Contact-angle boundary condition for multi-phase interface-capturing
+    simulations.  Used in conjuction with multiphaseMixture.
+
+SourceFiles
+    alphaContactAngleFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef alphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+#include "multiphaseMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class alphaContactAngleFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class alphaContactAngleFvPatchScalarField
+:
+    public zeroGradientFvPatchScalarField
+{
+public:
+
+    class interfaceThetaProps
+    {
+        //- Equilibrium contact angle
+        scalar theta0_;
+
+        //- Dynamic contact angle velocity scale
+        scalar uTheta_;
+
+        //- Limiting advancing contact angle
+        scalar thetaA_;
+
+        //- Limiting receeding contact angle
+        scalar thetaR_;
+
+
+    public:
+
+        // Constructors
+            interfaceThetaProps()
+            {}
+
+            interfaceThetaProps(Istream&);
+
+
+        // Member functions
+
+            //- Return the equilibrium contact angle theta0
+            scalar theta0(bool matched=true) const
+            {
+                if (matched) return theta0_;
+                else return 180.0 - theta0_;
+            }
+
+            //- Return the dynamic contact angle velocity scale
+            scalar uTheta() const
+            {
+                return uTheta_;
+            }
+
+            //- Return the limiting advancing contact angle
+            scalar thetaA(bool matched=true) const
+            {
+                if (matched) return thetaA_;
+                else return 180.0 - thetaA_;
+            }
+
+            //- Return the limiting receeding contact angle
+            scalar thetaR(bool matched=true) const
+            {
+                if (matched) return thetaR_;
+                else return 180.0 - thetaR_;
+            }
+
+
+        // IO functions
+
+            friend Istream& operator>>(Istream&, interfaceThetaProps&);
+            friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
+    };
+
+    typedef HashTable
+    <
+        interfaceThetaProps,
+        multiphaseMixture::interfacePair,
+        multiphaseMixture::interfacePair::hash
+    > thetaPropsTable;
+
+
+private:
+
+    // Private data
+
+        thetaPropsTable thetaProps_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("alphaContactAngle");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given alphaContactAngleFvPatchScalarField
+        //  onto a new patch
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return the contact angle properties
+        const thetaPropsTable& thetaProps() const
+        {
+            return thetaProps_;
+        }
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C

@@ -0,0 +1,772 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiphaseMixture.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "Time.H"
+#include "subCycle.H"
+#include "MULES.H"
+#include "surfaceInterpolate.H"
+#include "fvcGrad.H"
+#include "fvcSnGrad.H"
+#include "fvcDiv.H"
+#include "fvcFlux.H"
+
+// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
+
+const Foam::scalar Foam::multiphaseMixture::convertToRad =
+    Foam::constant::mathematical::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::multiphaseMixture::calcAlphas()
+{
+    scalar level = 0.0;
+    alphas_ == 0.0;
+
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        alphas_ += level*iter();
+        level += 1.0;
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::calcNu() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<volScalarField> tnu = iter()*iter().nu();
+    volScalarField& nu = tnu.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        nu += iter()*iter().nu();
+    }
+
+    return tnu;
+}
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixture::calcStf() const
+{
+    tmp<surfaceScalarField> tstf
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "stf",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "stf",
+                dimensionSet(1, -2, -2, 0, 0),
+                0.0
+            )
+        )
+    );
+
+    surfaceScalarField& stf = tstf.ref();
+
+    forAllConstIter(PtrDictionary<phase>, phases_, iter1)
+    {
+        const phase& alpha1 = iter1();
+
+        PtrDictionary<phase>::const_iterator iter2 = iter1;
+        ++iter2;
+
+        for (; iter2 != phases_.end(); ++iter2)
+        {
+            const phase& alpha2 = iter2();
+
+            sigmaTable::const_iterator sigma =
+                sigmas_.find(interfacePair(alpha1, alpha2));
+
+            if (sigma == sigmas_.end())
+            {
+                FatalErrorInFunction
+                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
+                    << " in list of sigma values"
+                    << exit(FatalError);
+            }
+
+            stf += dimensionedScalar("sigma", dimSigma_, sigma())
+               *fvc::interpolate(K(alpha1, alpha2))*
+                (
+                    fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
+                  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
+                );
+        }
+    }
+
+    return tstf;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiphaseMixture::multiphaseMixture
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi,
+    const volScalarField& voidfraction
+)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "transportProperties",
+            U.time().constant(),
+            U.db(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+
+    phases_(lookup("phases"), phase::iNew(U, phi)),
+
+    mesh_(U.mesh()),
+    U_(U),
+    phi_(phi),
+    voidfraction_(voidfraction),
+    rhoPhi_
+    (
+        IOobject
+        (
+            "rhoPhi",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("rhoPhi", dimMass/dimTime, 0.0)
+    ),
+    surfaceTensionForce_
+    (
+        IOobject
+        (
+            "surfaceTensionForce",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("surfaceTensionForce", dimensionSet(1, -2, -2, 0, 0), 0.0)
+    ),
+    alphas_
+    (
+        IOobject
+        (
+            "alphas",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("alphas", dimless, 0.0)
+    ),
+
+    nu_
+    (
+        IOobject
+        (
+            "nu",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        calcNu()
+    ),
+
+    sigmas_(lookup("sigmas")),
+    dimSigma_(1, 0, -2, 0, 0),
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mesh_.V()), 1.0/3.0)
+    )
+{
+    calcAlphas();
+    alphas_.write();
+    surfaceTensionForce_ = calcStf();
+
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::rho() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<volScalarField> trho = iter()*iter().rho();
+    volScalarField& rho = trho.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        rho += iter()*iter().rho();
+    }
+
+    return trho;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::multiphaseMixture::rho(const label patchi) const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
+    scalarField& rho = trho.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        rho += iter().boundaryField()[patchi]*iter().rho().value();
+    }
+
+    return trho;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::mu() const
+{
+    return rho()*nu();
+//    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+//    tmp<volScalarField> tmu = iter()*iter().rho()*iter().nu();
+//    volScalarField& mu = tmu.ref();
+
+//    for (++iter; iter != phases_.end(); ++iter)
+//    {
+//        mu += iter()*iter().rho()*iter().nu();
+//    }
+
+//    return tmu;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::multiphaseMixture::mu(const label patchi) const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<scalarField> tmu =
+        iter().boundaryField()[patchi]
+       *iter().rho().value()
+       *iter().nu(patchi);
+    scalarField& mu = tmu.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        mu +=
+            iter().boundaryField()[patchi]
+           *iter().rho().value()
+           *iter().nu(patchi);
+    }
+
+    return tmu;
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixture::muf() const
+{
+
+    return nuf()*fvc::interpolate(rho());
+//    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+//    tmp<surfaceScalarField> tmuf =
+//            fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
+//    surfaceScalarField& muf = tmuf.ref();
+
+//    for (++iter; iter != phases_.end(); ++iter)
+//    {
+//        muf +=
+//                fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
+//    }
+
+//    return tmuf;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::nu() const
+{
+    return nu_;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::multiphaseMixture::nu(const label patchi) const
+{
+    //return nu_.boundaryField()[patchi];
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<scalarField> tnu =
+        iter().boundaryField()[patchi]
+       *iter().nu(patchi);
+    scalarField& nu = tnu.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        nu +=
+            iter().boundaryField()[patchi]
+           *iter().nu(patchi);
+    }
+
+    return tnu;
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixture::nuf() const
+{
+    //return muf()/fvc::interpolate(rho());
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<surfaceScalarField> tnuf =
+        fvc::interpolate(iter())*fvc::interpolate(iter().nu());
+    surfaceScalarField& nuf = tnuf.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        nuf +=
+            fvc::interpolate(iter())*fvc::interpolate(iter().nu());
+    }
+
+    return tnuf;
+}
+
+void Foam::multiphaseMixture::solve()
+{
+    correct();
+
+    const Time& runTime = mesh_.time();
+
+    volScalarField& alpha = phases_.first();
+
+    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+    scalar cAlpha(readScalar(alphaControls.lookup("cAlpha")));
+
+    if (nAlphaSubCycles > 1)
+    {
+        surfaceScalarField rhoPhiSum
+        (
+            IOobject
+            (
+                "rhoPhiSum",
+                runTime.timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("0", rhoPhi_.dimensions(), 0)
+        );
+
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            FatalError << "Sub-cycling of the alpha equation not yet implemented!!" << abort(FatalError);
+            solveAlphas(cAlpha);
+            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
+        }
+
+        rhoPhi_ = rhoPhiSum;
+    }
+    else
+    {
+        solveAlphas(cAlpha);
+    }
+
+    // Update the mixture kinematic viscosity
+    nu_ = calcNu();
+    surfaceTensionForce_ = calcStf();
+}
+
+
+void Foam::multiphaseMixture::correct()
+{
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        iter().correct();
+    }
+}
+
+
+Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    /*
+    // Cell gradient of alpha
+    volVectorField gradAlpha =
+        alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
+
+    // Interpolated face-gradient of alpha
+    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
+    */
+
+    surfaceVectorField gradAlphaf
+    (
+        fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
+      - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
+    );
+
+    // Face unit interface normal
+    return gradAlphaf/(mag(gradAlphaf) + deltaN_);
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    // Face unit interface normal flux
+    return nHatfv(alpha1, alpha2) & mesh_.Sf();
+}
+
+
+// Correction for the boundary condition on the unit normal nHat on
+// walls to produce the correct contact angle.
+
+// The dynamic contact angle is calculated from the component of the
+// velocity on the direction of the interface, parallel to the wall.
+
+void Foam::multiphaseMixture::correctContactAngle
+(
+    const phase& alpha1,
+    const phase& alpha2,
+    surfaceVectorField::Boundary& nHatb
+) const
+{
+    const volScalarField::Boundary& gbf
+        = alpha1.boundaryField();
+
+    const fvBoundaryMesh& boundary = mesh_.boundary();
+
+    forAll(boundary, patchi)
+    {
+        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
+        {
+            const alphaContactAngleFvPatchScalarField& acap =
+                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
+
+            vectorField& nHatPatch = nHatb[patchi];
+
+            vectorField AfHatPatch
+            (
+                mesh_.Sf().boundaryField()[patchi]
+               /mesh_.magSf().boundaryField()[patchi]
+            );
+
+            alphaContactAngleFvPatchScalarField::thetaPropsTable::
+                const_iterator tp =
+                acap.thetaProps().find(interfacePair(alpha1, alpha2));
+
+            if (tp == acap.thetaProps().end())
+            {
+                FatalErrorInFunction
+                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
+                    << "\n    in table of theta properties for patch "
+                    << acap.patch().name()
+                    << exit(FatalError);
+            }
+
+            bool matched = (tp.key().first() == alpha1.name());
+
+            scalar theta0 = convertToRad*tp().theta0(matched);
+            scalarField theta(boundary[patchi].size(), theta0);
+
+            scalar uTheta = tp().uTheta();
+
+            // Calculate the dynamic contact angle if required
+            if (uTheta > SMALL)
+            {
+                scalar thetaA = convertToRad*tp().thetaA(matched);
+                scalar thetaR = convertToRad*tp().thetaR(matched);
+
+                // Calculated the component of the velocity parallel to the wall
+                vectorField Uwall
+                (
+                    U_.boundaryField()[patchi].patchInternalField()
+                  - U_.boundaryField()[patchi]
+                );
+                Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
+
+                // Find the direction of the interface parallel to the wall
+                vectorField nWall
+                (
+                    nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
+                );
+
+                // Normalise nWall
+                nWall /= (mag(nWall) + SMALL);
+
+                // Calculate Uwall resolved normal to the interface parallel to
+                // the interface
+                scalarField uwall(nWall & Uwall);
+
+                theta += (thetaA - thetaR)*tanh(uwall/uTheta);
+            }
+
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12(nHatPatch & AfHatPatch);
+
+            scalarField b1(cos(theta));
+
+            scalarField b2(nHatPatch.size());
+
+            forAll(b2, facei)
+            {
+                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
+            }
+
+            scalarField det(1.0 - a12*a12);
+
+            scalarField a((b1 - a12*b2)/det);
+            scalarField b((b2 - a12*b1)/det);
+
+            nHatPatch = a*AfHatPatch + b*nHatPatch;
+
+            nHatPatch /= (mag(nHatPatch) + deltaN_.value());
+        }
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
+(
+    const phase& alpha1,
+    const phase& alpha2
+) const
+{
+    tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
+
+    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
+
+    // Simple expression for curvature
+    return -fvc::div(tnHatfv & mesh_.Sf());
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::nearInterface() const
+{
+    tmp<volScalarField> tnearInt
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "nearInterface",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("nearInterface", dimless, 0.0)
+        )
+    );
+
+    forAllConstIter(PtrDictionary<phase>, phases_, iter)
+    {
+        tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
+    }
+
+    return tnearInt;
+}
+
+
+void Foam::multiphaseMixture::solveAlphas
+(
+    const scalar cAlpha
+)
+{
+    static label nSolves=-1;
+    nSolves++;
+
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phic(mag(phi_/mesh_.magSf()));
+    phic = min(cAlpha*phic, max(phic));
+
+    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
+    int phasei = 0;
+
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        phase& alpha = iter();
+
+        alphaPhiCorrs.set
+        (
+            phasei,
+            new surfaceScalarField
+            (
+                "phi" + alpha.name() + "Corr",
+                fvc::flux
+                (
+                    phi_,
+                    alpha,
+                    alphaScheme
+                )
+            )
+        );
+
+        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
+
+        forAllIter(PtrDictionary<phase>, phases_, iter2)
+        {
+            phase& alpha2 = iter2();
+
+            if (&alpha2 == &alpha) continue;
+
+            surfaceScalarField phir(phic*nHatf(alpha, alpha2));
+
+            alphaPhiCorr += fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha,
+                alpharScheme
+            );
+        }
+
+        MULES::limit
+        (
+            1.0/mesh_.time().deltaT().value(),
+            voidfraction_,
+            alpha,
+            phi_,
+            alphaPhiCorr,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            true
+        );
+
+        phasei++;
+    }
+
+    MULES::limitSum(alphaPhiCorrs);
+
+    rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
+
+    volScalarField sumAlpha
+    (
+        IOobject
+        (
+            "sumAlpha",
+            mesh_.time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensionedScalar("sumAlpha", dimless, 0)
+    );
+
+    phasei = 0;
+
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        phase& alpha = iter();
+
+        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
+
+        MULES::explicitSolve
+        (
+            voidfraction_,
+            alpha,
+            alphaPhi,
+            zeroField(),
+            zeroField()
+        );
+
+        rhoPhi_ += alphaPhi*alpha.rho();
+
+        Info<< alpha.name() << " volume fraction, min, max = "
+            << alpha.weightedAverage(mesh_.V()).value()
+            << ' ' << min(alpha).value()
+            << ' ' << max(alpha).value()
+            << endl;
+
+        sumAlpha += alpha;
+
+        phasei++;
+    }
+
+    Info<< "Phase-sum volume fraction, min, max = "
+        << sumAlpha.weightedAverage(mesh_.V()).value()
+        << ' ' << min(sumAlpha).value()
+        << ' ' << max(sumAlpha).value()
+        << endl;
+
+    calcAlphas();
+}
+
+
+bool Foam::multiphaseMixture::read()
+{
+    if (transportModel::read())
+    {
+        bool readOK = true;
+
+        PtrList<entry> phaseData(lookup("phases"));
+        label phasei = 0;
+
+        forAllIter(PtrDictionary<phase>, phases_, iter)
+        {
+            readOK &= iter().read(phaseData[phasei++].dict());
+        }
+
+        lookup("sigmas") >> sigmas_;
+
+        return readOK;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.H

@@ -0,0 +1,284 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Class
+    multiphaseMixture
+
+Description
+    This class is based on the OpenFOAM(R) Foam::multiphaseMixture class,
+    which is an incompressible multi-phase mixture with built in solution
+    for the phase fractions with interface compression for interface-capturing.
+    It has been extended to include the void fraction in the volume fraction
+    transport equations.
+
+    Derived from transportModel so that it can be unsed in conjunction with
+    the incompressible turbulence models.
+
+    Surface tension and contact-angle is handled for the interface
+    between each phase-pair.
+
+SourceFiles
+    multiphaseMixture.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMixture_H
+#define multiphaseMixture_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "IOdictionary.H"
+#include "phase.H"
+#include "PtrDictionary.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class multiphaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMixture
+:
+    public IOdictionary,
+    public transportModel
+{
+public:
+
+    class interfacePair
+    :
+        public Pair<word>
+    {
+    public:
+
+        class hash
+        :
+            public Hash<interfacePair>
+        {
+        public:
+
+            hash()
+            {}
+
+            label operator()(const interfacePair& key) const
+            {
+                return word::hash()(key.first()) + word::hash()(key.second());
+            }
+        };
+
+
+        // Constructors
+
+            interfacePair()
+            {}
+
+            interfacePair(const word& alpha1Name, const word& alpha2Name)
+            :
+                Pair<word>(alpha1Name, alpha2Name)
+            {}
+
+            interfacePair(const phase& alpha1, const phase& alpha2)
+            :
+                Pair<word>(alpha1.name(), alpha2.name())
+            {}
+
+
+        // Friend Operators
+
+            friend bool operator==
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return
+                (
+                    ((a.first() == b.first()) && (a.second() == b.second()))
+                 || ((a.first() == b.second()) && (a.second() == b.first()))
+                );
+            }
+
+            friend bool operator!=
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return (!(a == b));
+            }
+    };
+
+
+private:
+
+    // Private data
+
+        //- Dictionary of phases
+        PtrDictionary<phase> phases_;
+
+        const fvMesh& mesh_;
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+        const volScalarField& voidfraction_;
+        surfaceScalarField rhoPhi_;
+        surfaceScalarField surfaceTensionForce_;
+        volScalarField alphas_;
+
+        volScalarField nu_;
+
+        typedef HashTable<scalar, interfacePair, interfacePair::hash>
+            sigmaTable;
+
+        sigmaTable sigmas_;
+        dimensionSet dimSigma_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        //- Conversion factor for degrees into radians
+        static const scalar convertToRad;
+
+
+    // Private member functions
+
+        void calcAlphas();
+
+        tmp<volScalarField> calcNu() const;
+
+        void solveAlphas(const scalar cAlpha);
+
+        tmp<surfaceVectorField> nHatfv
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<surfaceScalarField> nHatf
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        void correctContactAngle
+        (
+            const phase& alpha1,
+            const phase& alpha2,
+            surfaceVectorField::Boundary& nHatb
+        ) const;
+
+        tmp<volScalarField> K(const phase& alpha1, const phase& alpha2) const;
+        tmp<surfaceScalarField> calcStf() const;
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        multiphaseMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            const volScalarField& voidfraction
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMixture()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phases
+        const PtrDictionary<phase>& phases() const
+        {
+            return phases_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the volumetric flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        const surfaceScalarField& rhoPhi() const
+        {
+            return rhoPhi_;
+        }
+
+        //- Return the mixture density
+        tmp<volScalarField> rho() const;
+
+        //- Return the mixture density for patch
+        tmp<scalarField> rho(const label patchi) const;
+
+        //- Return the dynamic laminar viscosity
+        tmp<volScalarField> mu() const;
+
+        //- Return the dynamic laminar viscosity for patch
+        tmp<scalarField> mu(const label patchi) const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> muf() const;
+
+        //- Return the kinematic laminar viscosity
+        tmp<volScalarField> nu() const;
+
+        //- Return the laminar viscosity for patch
+        tmp<scalarField> nu(const label patchi) const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> nuf() const;
+
+        tmp<surfaceScalarField> surfaceTensionForce() const
+        {
+            return surfaceTensionForce_;
+        }
+
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        //- Solve for the mixture phase-fractions
+        void solve();
+
+        //- Correct the mixture properties
+        void correct();
+
+        //- Read base transportProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/phase/phase.C

@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phase.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phase::phase
+(
+    const word& phaseName,
+    const dictionary& phaseDict,
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    volScalarField
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phaseName),
+            U.mesh().time().timeName(),
+            U.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U.mesh()
+    ),
+    name_(phaseName),
+    phaseDict_(phaseDict),
+    nuModel_
+    (
+        viscosityModel::New
+        (
+            IOobject::groupName("nu", phaseName),
+            phaseDict_,
+            U,
+            phi
+        )
+    ),
+    rho_("rho", dimDensity, phaseDict_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::phase> Foam::phase::clone() const
+{
+    NotImplemented;
+    return autoPtr<phase>(NULL);
+}
+
+
+void Foam::phase::correct()
+{
+    nuModel_->correct();
+}
+
+
+bool Foam::phase::read(const dictionary& phaseDict)
+{
+    phaseDict_ = phaseDict;
+
+    if (nuModel_->read(phaseDict_))
+    {
+        phaseDict_.lookup("rho") >> rho_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/phase/phase.H

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phase
+
+Description
+    Single incompressible phase derived from the phase-fraction.
+    Used as part of the multiPhaseMixture for interface-capturing multi-phase
+    simulations.
+
+SourceFiles
+    phase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phase_H
+#define phase_H
+
+#include "volFields.H"
+#include "dictionaryEntry.H"
+#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class phase Declaration
+\*---------------------------------------------------------------------------*/
+
+class phase
+:
+    public volScalarField
+{
+    // Private data
+
+        word name_;
+        dictionary phaseDict_;
+        autoPtr<viscosityModel> nuModel_;
+        dimensionedScalar rho_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        phase
+        (
+            const word& name,
+            const dictionary& phaseDict,
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+        //- Return clone
+        autoPtr<phase> clone() const;
+
+        //- Return a pointer to a new phase created on freestore
+        //  from Istream
+        class iNew
+        {
+            const volVectorField& U_;
+            const surfaceScalarField& phi_;
+
+        public:
+
+            iNew
+            (
+                const volVectorField& U,
+                const surfaceScalarField& phi
+            )
+            :
+                U_(U),
+                phi_(phi)
+            {}
+
+            autoPtr<phase> operator()(Istream& is) const
+            {
+                dictionaryEntry ent(dictionary::null, is);
+                return autoPtr<phase>(new phase(ent.keyword(), ent, U_, phi_));
+            }
+        };
+
+
+    // Member Functions
+
+        const word& name() const
+        {
+            return name_;
+        }
+
+        const word& keyword() const
+        {
+            return name();
+        }
+
+        //- Return const-access to phase1 viscosityModel
+        const viscosityModel& nuModel() const
+        {
+            return nuModel_();
+        }
+
+        //- Return the kinematic laminar viscosity
+        tmp<volScalarField> nu() const
+        {
+            return nuModel_->nu();
+        }
+
+        //- Return the laminar viscosity for patch
+        tmp<scalarField> nu(const label patchi) const
+        {
+            return nuModel_->nu(patchi);
+        }
+
+        //- Return const-access to phase1 density
+        const dimensionedScalar& rho() const
+        {
+            return rho_;
+        }
+
+        //- Correct the phase properties
+        void correct();
+
+        //-Inherit read from volScalarField
+        using volScalarField::read;
+
+        //- Read base transportProperties dictionary
+        bool read(const dictionary& phaseDict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/pEqn.H

@@ -0,0 +1,73 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUepsf("rAUepsf", fvc::interpolate(rAU*voidfraction));
+    surfaceScalarField rAUepsSqf("rAUepsSqf", fvc::interpolate(rAU*voidfraction*voidfraction));
+    volVectorField Ueps("Ueps", U * voidfraction);
+
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA*voidfraction)
+      + fvc::interpolate(voidfraction*rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+
+    adjustPhi(phiHbyA, U, p_rgh);
+
+    if (modelType == "A")
+        rAUepsf = rAUepsSqf;
+
+    surfaceScalarField phig (-ghf*fvc::snGrad(rho)*rAUepsf*mesh.magSf());
+
+    surfaceScalarField phiSt (mixture.surfaceTensionForce()*rAUepsSqf*mesh.magSf());
+
+    surfaceScalarField phiS (fvc::flux(voidfraction*Us*Ksl*rAU));
+
+    phiHbyA += phig + phiSt + phiS;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, Ueps, phiHbyA, rAUepsf);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUepsf, p_rgh) == particleCloud.ddtVoidfraction() + fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            if (modelType == "A")
+                U = HbyA + voidfraction*rAU*fvc::reconstruct((phig-p_rghEqn.flux()+phiSt)/rAUepsf) + rAU*Us*Ksl;
+            else
+                U = HbyA + rAU*fvc::reconstruct((phig-p_rghEqn.flux()+phiSt)/rAUepsf) + rAU*Us*Ksl;
+
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+}

applications/solvers/cfdemSolverPiso/UEqn.H

@@ -15,10 +15,10 @@ fvOptions.constrain(UEqn);
 if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
 {
     solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-    fvOptions.correct(U);		    
+    fvOptions.correct(U);
 }
 else if (piso.momentumPredictor())
 {
     solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
     fvOptions.correct(U);
-}
+}

applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C

@@ -86,12 +86,12 @@ int main(int argc, char *argv[])
         Ksl = particleCloud.momCoupleM(0).impMomSource();
         Ksl.correctBoundaryConditions();
 
-       //Force Checks
-       vector fTotal(0,0,0);
-       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
-       reduce(fImpTotal, sumOp<vector>());
-       Info << "TotalForceExp: " << fTotal << endl;
-       Info << "TotalForceImp: " << fImpTotal << endl;
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
 
         #include "solverDebugInfo.H"
         particleCloud.clockM().stop("Coupling");

applications/solvers/cfdemSolverPiso/createFields.H

@@ -96,17 +96,17 @@
 #define createPhi_H
 Info<< "Reading/calculating face flux field phi\n" << endl;
 surfaceScalarField phi
- (
-     IOobject
-     (
+(
+    IOobject
+    (
         "phi",
         runTime.timeName(),
          mesh,
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
-     ),
-     linearInterpolate(U*voidfraction) & mesh.Sf()
- );
+    ),
+    linearInterpolate(U*voidfraction) & mesh.Sf()
+);
 #endif
 
 
@@ -123,4 +123,4 @@ surfaceScalarField phi
         incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
 
-#include "createMRF.H"
+#include "createMRF.H"

applications/solvers/cfdemSolverPiso/pEqn.H

@@ -31,12 +31,12 @@ constrainPressure(p, Uvoidfraction, phiHbyA, rAUvoidfraction, MRF);
 while (piso.correctNonOrthogonal())
 {
     // Pressure corrector
-  
+
     fvScalarMatrix pEqn
     (
         fvm::laplacian(rAUvoidfraction, p) == fvc::div(phi) + particleCloud.ddtVoidfraction()
     );
-    
+
     pEqn.setReference(pRefCell, pRefValue);
 
     pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
@@ -55,4 +55,4 @@ else
     U = HbyA - voidfraction*rAU*fvc::grad(p) + Ksl/rho*Us*rAU;
 
 U.correctBoundaryConditions();
-fvOptions.correct(U);
+fvOptions.correct(U);

applications/solvers/cfdemSolverPisoScalar/TEqn.H

@@ -1,4 +1,4 @@
- // get scalar source from DEM        
+ // get scalar source from DEM
         particleCloud.forceM(1).manipulateScalarField(Tsource);
         Tsource.correctBoundaryConditions();
 
@@ -12,4 +12,4 @@
            Tsource
         );
         TEqn.relax();
-        TEqn.solve();
+        TEqn.solve();

applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C

@@ -81,23 +81,23 @@ int main(int argc, char *argv[])
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
         }
-    
+
         Info << "update Ksl.internalField()" << endl;
         Ksl = particleCloud.momCoupleM(0).impMomSource();
         Ksl.correctBoundaryConditions();
 
-       //Force Checks
-       vector fTotal(0,0,0);
-       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
-       reduce(fImpTotal, sumOp<vector>());
-       Info << "TotalForceExp: " << fTotal << endl;
-       Info << "TotalForceImp: " << fImpTotal << endl;
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
 
         #include "solverDebugInfo.H"
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
-	
+
         #include "TEqn.H"
 
         if(particleCloud.solveFlow())

applications/solvers/cfdemSolverPisoScalar/createFields.H

@@ -146,17 +146,17 @@
 #define createPhi_H
 Info<< "Reading/calculating face flux field phi\n" << endl;
 surfaceScalarField phi
- (
-     IOobject
-     (
+(
+    IOobject
+    (
         "phi",
         runTime.timeName(),
          mesh,
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
-     ),
-     linearInterpolate(U*voidfraction) & mesh.Sf()
- );
+    ),
+    linearInterpolate(U*voidfraction) & mesh.Sf()
+);
 #endif
 
 
@@ -173,4 +173,4 @@ surfaceScalarField phi
         incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
 
-#include "createMRF.H"
+#include "createMRF.H"

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -23,6 +23,9 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+    // correct source for the thermodynamic reference temperature
+    dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
+    Qsource += QCoeff*Tref;
 
     fvScalarMatrix EEqn
     (

applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C

@@ -69,8 +69,8 @@ int main(int argc, char *argv[])
     #include "checkModelType.H"
 
     turbulence->validate();
-  //        #include "compressibleCourantNo.H"
-  //  #include "setInitialDeltaT.H"
+    //#include "compressibleCourantNo.H"
+    //#include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/cfdemSolverRhoPimple/pEqn.H

@@ -32,7 +32,7 @@ else
     //      + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
-    
+
     // flux without pressure gradient contribution
     phi = phiHbyA + phiUs;
 

applications/solvers/cfdemSolverRhoPimple/rhoEqn.H

@@ -14,4 +14,4 @@
     fvOptions.correct(rho);
 }
 
-// ************************************************************************* //
+// ************************************************************************* //

applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H

@@ -9,6 +9,10 @@ particleCloud.energyCoefficients(QCoeff);
 thCond=particleCloud.thermCondM().thermCond();
 Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+// correct source for the thermodynamic reference temperature
+dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
+Qsource += QCoeff*Tref;
+
 fvScalarMatrix EEqn
 (
     fvm::ddt(rhoeps, he) + fvm::div(phi, he)

applications/solvers/cfdemSolverRhoSimple/EEqn.H

@@ -6,13 +6,13 @@
     particleCloud.energyContributions(Qsource);
     particleCloud.energyCoefficients(QCoeff);
 
-    //thDiff=particleCloud.thermCondM().thermDiff();   
-    thCond=particleCloud.thermCondM().thermCond(); 
+    //thDiff=particleCloud.thermCondM().thermDiff();
+    thCond=particleCloud.thermCondM().thermCond();
 
-    addSource = 
+    addSource =
         (
             he.name() == "e"
-          ? 
+          ?
             fvc::div(phi, K) +
             fvc::div
             (
@@ -25,6 +25,9 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+    // correct source for the thermodynamic reference temperature
+    dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
+    Qsource += QCoeff*Tref;
 
     fvScalarMatrix EEqn
     (

applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C

@@ -17,14 +17,12 @@ License
     Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
 
 Application
-    cfdemSolverRhoPimple
+    cfdemSolverRhoSimple
 
 Description
-    Transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
-    algorithm.
-    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
-    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
-    where additional functionality for CFD-DEM coupling is added.
+    Steady-state solver for turbulent flow of compressible fluids based on
+    rhoSimpleFoam where functionality for CFD-DEM coupling has been added.
+
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"

applications/solvers/cfdemSolverRhoSimple/pEqn.H

@@ -27,7 +27,7 @@ else
             fvc::flux(rhoeps*HbyA)
         )
     );
-    
+
     // flux without pressure gradient contribution
     phi = phiHbyA + phiUs;
 
@@ -78,4 +78,4 @@ else
 
 U.correctBoundaryConditions();
 fvOptions.correct(U);
-K = 0.5*magSqr(U);
+K = 0.5*magSqr(U);

applications/utilities/cfdemPostproc/cfdemPostproc.C

@@ -84,7 +84,7 @@ int main(int argc, char *argv[])
     particleCloud.dataExchangeM().allocateArray(particleV_,0.,1);
     particleCloud.get_cellIDs(cellIDs_);  // get ref to cellIDs
     //particleCloud.dataExchangeM().allocateArray(cellIDs_,0.,1);
-    
+
 
     while (runTime.loop())
     {

applications/utilities/writeUfluid/writeUfluid.C

@@ -28,7 +28,7 @@ Application
     writeUfluidwriteUfluid
 
 Description
-    Writes the the cell center fluid velocity to particles in the lagrangian 
+    Writes the the cell center fluid velocity to particles in the lagrangian
     time directory.
 \*---------------------------------------------------------------------------*/
 
@@ -76,13 +76,13 @@ int nParticle=0;
     {
         volVectorField U(UHeader,mesh);
         passiveParticleCloud myCloud(mesh, cloudName);
-		myCloud.write();
+        myCloud.write();
         nParticle = myCloud.size();
-	    IOField<vector> Ufluid(myCloud.fieldIOobject("Ufluid",IOobject::NO_READ),nParticle);
+        IOField<vector> Ufluid(myCloud.fieldIOobject("Ufluid",IOobject::NO_READ),nParticle);
         label i = 0;
         forAllConstIter(passiveParticleCloud, myCloud, iter)
         {
-		    Ufluid[i]=U[iter().cell()];
+            Ufluid[i]=U[iter().cell()];
             i++;
         }
         Ufluid.write();

doc/CFDEMcoupling_install.txt

@@ -189,7 +189,7 @@ cfdemCompLIG :pre
 
 If the compilation fails with a message like
 
-No rule to make target `/usr/lib/libpython2.7.so' :pre
+No rule to make target '/usr/lib/libpython2.7.so' :pre
 
 you probably need to create a symbolic link to the library in question.
 

doc/CFDEMcoupling_models.txt

@@ -38,7 +38,7 @@ models used for chemical reaction calculations.
 "diffusionCoefficients"_chemistryModel_diffusionCoefficients.html,
 "massTransferCoeff"_chemistryModel_massTransferCoeff.html,
 "off"_chemistryModel_noChemistry.html,
-reactantPerParticle,
+"reactantPerParticle"_chemistryModel_reactantPerParticle.html,
 "species"_chemistryModel_species.html :tb(c=2,ea=c)
 
 
@@ -60,7 +60,8 @@ that performs the data exchange between the DEM code and the CFD code.
 "oneWayVTK"_dataExchangeModel_oneWayVTK.html,
 "twoWayFiles"_dataExchangeModel_twoWayFiles.html,
 "twoWayMPI"_dataExchangeModel_twoWayMPI.html,
-"twoWayMany2Many"_dataExchangeModel_twoWayMany2Many.html :tb(c=2,ea=c)
+"twoWayMany2Many"_dataExchangeModel_twoWayMany2Many.html,
+"twoWayOne2One"_dataExchangeModel_twoWayOne2One.html :tb(c=2,ea=c)
 
 
 6.6 Energy models :h4
@@ -94,11 +95,13 @@ Fines,
 "fieldStore"_forceModel_fieldStore.html,
 "fieldTimeAverage"_forceModel_fieldTimeAverage.html,
 "gradPForce"_forceModel_gradPForce.html,
+"gradPForceSmooth"_forceModel_gradPForceSmooth.html,
 granKineticEnergy,
 "interface"_forceModel_interface.html,
 "noDrag"_forceModel_noDrag.html,
 "particleCellVolume"_forceModel_particleCellVolume.html,
 "pdCorrelation"_forceModel_pdCorrelation.html,
+"surfaceTensionForce"_forceModel_surfaceTensionForce.html,
 "virtualMassForce"_forceModel_virtualMassForce.html,
 "viscForce"_forceModel_viscForce.html,
 "volWeightedAverage"_forceModel_volWeightedAverage.html :tb(c=2,ea=c)
@@ -188,7 +191,8 @@ The "smoothingModel"_smoothingModel.html keyword entry specifies the model for
 smoothing the exchange fields.
 
 "constDiffSmoothing"_smoothingModel_constDiffSmoothing.html,
-"off"_smoothingModel_noSmoothing.html :tb(c=2,ea=c)
+"off"_smoothingModel_noSmoothing.html,
+"temporalSmoothing"_smoothingModel_temporalSmoothing.html :tb(c=2,ea=c)
 
 
 6.16 Thermal conductivity models :h4

doc/CFDEMcoupling_solvers.txt

@@ -10,9 +10,10 @@
 This section lists all CFDEMcoupling solvers alphabetically.
 
 "cfdemSolverIB"_cfdemSolverIB.html,
+"cfdemSolverMultiphase"_cfdemSolverMultiphase.html,
 "cfdemSolverPiso"_cfdemSolverPiso.html,
 "cfdemSolverPisoScalar"_cfdemSolverPisoScalar.html,
-cfdemSolverRhoPimple,
-cfdemSolverRhoPimpleChem,
-cfdemSolverRhoSimple :tb(c=2,ea=c)
+"cfdemSolverRhoPimple"_cfdemSolverRhoPimple.html,
+"cfdemSolverRhoPimpleChem"_cfdemSolverRhoPimpleChem.html,
+"cfdemSolverRhoSimple"_cfdemSolverRhoSimple.html :tb(c=2,ea=c)
 

doc/cfdemSolverMultiphase.txt

@@ -0,0 +1,62 @@
+<!-- HTML_ONLY -->
+<HEAD>
+  <META CHARSET="utf-8">
+</HEAD>
+<!-- END_HTML_ONLY -->
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+cfdemSolverMultiphase command :h3
+
+[Description:]
+
+<!-- HTML_ONLY -->
+"cfdemSolverMultiphase" is a coupled CFD-DEM solver using the CFDEMcoupling framework. Based on the OpenFOAM solver multiphaseInterFoam&reg;(*) it has functionality to simulate several fluids using the Volume of Fluid approach, coupled with the DEM code LIGGGHTS for solid particles.
+
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+"cfdemSolverMultiphase" is a coupled CFD-DEM solver using the CFDEMcoupling framework. Based on the OpenFOAM solver multiphaseInterFoam\ |reg|\ (*) it has functionality to simulate several fluids using the Volume of Fluid approach, coupled with the DEM code LIGGGHTS for solid particles.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+
+For more details, see "Vångö et al. (2018)"_#Vångö2018.
+
+:line
+
+:link(Vångö2018)
+[(Vångö2018)] M. Vångö, S. Pirker, T. Lichtenegger. (2018):
+"Unresolved CFD-DEM modeling of multiphase flow in densely packed particle beds",
+Applied Mathematical Modelling
+
+:line
+
+<!-- HTML_ONLY -->
+NOTE:
+(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com, and owner of the
+OPENFOAM&reg; and OpenCFD&reg; trade marks.
+OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+.. note::
+
+   (*) This offering is not approved or endorsed by OpenCFD Limited, producer
+   and distributor of the OpenFOAM software via www.openfoam.com, and owner of
+   the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
+   OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
+   distributor of the OpenFOAM software via www.openfoam.com.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->

doc/cfdemSolverRhoPimple.txt

@@ -0,0 +1,59 @@
+<!-- HTML_ONLY -->
+<HEAD>
+  <META CHARSET="utf-8">
+</HEAD>
+<!-- END_HTML_ONLY -->
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+cfdemSolverRhoPimple command :h3
+
+[Description:]
+
+<!-- HTML_ONLY -->
+"cfdemSolverRhoPimple" is a coupled CFD-DEM solver using the CFDEMcoupling
+framework. Based on the OpenFOAM&reg;(*) solver rhoPimpleFoam, this is a
+transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
+algorithm, coupled with the DEM code LIGGGHTS for solid particles.
+
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+"cfdemSolverRhoPimple" is a coupled CFD-DEM solver using the CFDEMcoupling
+framework. Based on the OpenFOAM\ |reg|\ (*) solver rhoPimpleFoam, this is a
+transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
+algorithm, coupled with the DEM code LIGGGHTS for solid particles.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+
+:line
+
+<!-- HTML_ONLY -->
+NOTE:
+(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com, and owner of the
+OPENFOAM&reg; and OpenCFD&reg; trade marks.
+OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+.. note::
+
+   (*) This offering is not approved or endorsed by OpenCFD Limited, producer
+   and distributor of the OpenFOAM software via www.openfoam.com, and owner of
+   the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
+   OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
+   distributor of the OpenFOAM software via www.openfoam.com.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->

doc/cfdemSolverRhoPimpleChem.txt

@@ -0,0 +1,63 @@
+<!-- HTML_ONLY -->
+<HEAD>
+  <META CHARSET="utf-8">
+</HEAD>
+<!-- END_HTML_ONLY -->
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+cfdemSolverRhoPimpleChem command :h3
+
+[Description:]
+
+<!-- HTML_ONLY -->
+"cfdemSolverRhoPimpleChem" is a coupled CFD-DEM solver using the CFDEMcoupling
+framework. Based on the OpenFOAM&reg;(*) solver rhoPimpleFoam, this is a
+transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
+algorithm, coupled with the DEM code LIGGGHTS for solid particles.
+Compared to cfdemSolverRhoPimple this solver adds functionality for chemical
+reactions.
+
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+"cfdemSolverRhoPimpleChem" is a coupled CFD-DEM solver using the CFDEMcoupling
+framework. Based on the OpenFOAM\ |reg|\ (*) solver rhoPimpleFoam, this is a
+transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
+algorithm, coupled with the DEM code LIGGGHTS for solid particles.
+Compared to cfdemSolverRhoPimple this solver adds functionality for chemical
+reactions.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+
+:line
+
+<!-- HTML_ONLY -->
+NOTE:
+(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com, and owner of the
+OPENFOAM&reg; and OpenCFD&reg; trade marks.
+OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+.. note::
+
+   (*) This offering is not approved or endorsed by OpenCFD Limited, producer
+   and distributor of the OpenFOAM software via www.openfoam.com, and owner of
+   the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
+   OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
+   distributor of the OpenFOAM software via www.openfoam.com.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->

doc/cfdemSolverRhoSimple.txt

@@ -0,0 +1,59 @@
+<!-- HTML_ONLY -->
+<HEAD>
+  <META CHARSET="utf-8">
+</HEAD>
+<!-- END_HTML_ONLY -->
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+cfdemSolverRhoSimple command :h3
+
+[Description:]
+
+<!-- HTML_ONLY -->
+"cfdemSolverRhoSimple" is a coupled CFD-DEM solver using the CFDEMcoupling
+framework. Based on the OpenFOAM&reg;(*) solver rhoSimpleFoam, this is a
+steady-state solver for turbulent flow of compressible fluids coupled with the
+DEM code LIGGGHTS for solid particles.
+
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+"cfdemSolverRhoSimple" is a coupled CFD-DEM solver using the CFDEMcoupling
+framework. Based on the OpenFOAM\ |reg|\ (*) solver rhoSimpleFoam, this is a
+steady-state solver for turbulent flow of compressible fluids coupled with the
+DEM code LIGGGHTS for solid particles.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+
+:line
+
+<!-- HTML_ONLY -->
+NOTE:
+(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com, and owner of the
+OPENFOAM&reg; and OpenCFD&reg; trade marks.
+OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+.. note::
+
+   (*) This offering is not approved or endorsed by OpenCFD Limited, producer
+   and distributor of the OpenFOAM software via www.openfoam.com, and owner of
+   the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
+   OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
+   distributor of the OpenFOAM software via www.openfoam.com.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->

doc/chemistryModel_diffusionCoefficients.txt

@@ -23,7 +23,7 @@ diffusionCoefficientsProps
     diffusantGasNames   ( speciesNames );
 \} :pre
 
-{switch1} = (optional, normally off) flag to give information :ulb,l
+{switch1} = (optional, default false) flag to output verbose information :ulb,l
 {ChemistryFile} = path to file, where the reacting species are listed :l
 {diffusantGasNames} = list of gas field names that are the reactant gases :l
 :ule

doc/chemistryModel_massTransferCoeff.txt

@@ -21,7 +21,7 @@ massTransferCoeffProps
     verbose         switch1;
 \} :pre
 
-{switch1} = (optional, normally off) flag to give information :l
+{switch1} = (optional, default false) flag to output verbose information :l
 :ule
 
 [Examples:]

doc/chemistryModel_reactantPerParticle.txt

@@ -0,0 +1,54 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+chemistryModel reactantPerParticle command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+chemistryModels
+(
+    reactantPerParticle
+);
+reactantPerParticleProps
+\{
+    voidfractionFieldName   "voidfraction";
+    Nevery                  number1;
+\} :pre
+
+{voidfraction} = (optional, default "voidfraction") name of the finite volume void fraction field :l
+{number1} = (optional, default 1) number to adjust execution interval :l
+:ule
+
+[Examples:]
+
+chemistryModels
+(
+    reactantPerParticle
+);
+reactantPerParticleProps
+\{
+    voidfractionFieldName   "voidfraction";
+    Nevery                  1;
+\} :pre
+
+[Description:]
+
+The chemistry model performs the calculation of chemical reactional effects
+acting on each DEM particle. The reactantPerParticle model is the model to
+communicate the available reactant per particle.
+
+[Restrictions:]
+
+none
+
+[Related commands:]
+
+"chemistryModel"_chemistryModel.html
+

doc/chemistryModel_species.txt

@@ -26,16 +26,18 @@ speciesProps
     partTempName            "partTemp";
     partRhoName             "partRho";
     verbose                 switch1;
+    Nevery                  number1;
 \} :pre
 
 {ChemistryFile} = path to file, where the reacting species are listed  :ulb,l
-{T} = name of the finite volume temperature field, it is already added in default and doesn't need to be specified if name is the same :l
-{rho} = name of the finite volume density field, it is already added in default and doesn't need to be specified if name is the same :l
-{voidfraction} = name of the finite volume void fraction field, it is already added in default and doesn't need to be specified if name is the same :l
-{molarConc} = name of the finite volume molar concentration field, it is already added in default and doesn't need to be specified if name is the same :l
-{partTemp} = name of the finite volume cell averaged particle temperature field, it is already added in default and doesn't need to be specified if name is the same :l
-{partRho} = name of the finite volume cell averaged density temperature field, it is already added in default and doesn't need to be specified if name is the same :l
-{switch1} = (optional, normally off) flag to give information :l
+{T} = (optional, default "T") name of the finite volume temperature field :l
+{rho} = (optional, default "rho") name of the finite volume density field :l
+{voidfraction} = (optional, default "voidfraction") name of the finite volume void fraction field :l
+{molarConc} = (optional, default "molarConc") name of the finite volume molar concentration field :l
+{partTemp} = (optional, default "partTemp") name of the finite volume cell averaged particle temperature field :l
+{partRho} = (optional, default "partRho") name of the finite volume cell averaged density temperature field :l
+{switch1} = (optional, default false) flag to output verbose information :l
+{number1} = (optional, default 1) number to adjust execution interval :l
 :ule
 
 [Examples:]

doc/dataExchangeModel_twoWayOne2One.txt

@@ -0,0 +1,47 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+dataExchangeModel_twoWayOne2One command :h3
+
+[Syntax:]
+
+Defined in couplingProperties dictionary.
+
+dataExchangeModel twoWayOne2One;
+twoWayOne2OneProps
+\{
+    liggghtsPath     "path";
+    useStaticProcMap switch1;
+\}; :pre
+
+{path}   = path to the DEM simulation input file :ulb,l
+{switch1} = (optional, default no) switch to determine if the map is built once (yes) or every coupling step (no) :l
+:ule
+
+[Examples:]
+
+dataExchangeModel twoWayOne2One;
+twoWayOne2OneProps
+\{
+    liggghtsPath     "../DEM/in.liggghts_init";
+    useStaticProcMap yes;
+\} :pre
+
+[Description:]
+
+The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayOne2One model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using a more sophisticated mapping scheme than twoWayMPI / all2all and scales much better for large systems and many cores. The DEM run is executed by the coupling model, via a liggghtsCommandModel object. Only use staticProcMap yes if no load balancing is employed.
+
+[Restrictions:]
+
+Must be used in combination with the engineSearchMany2Many locate model! Use the "one2one" cfd datacoupling option in fix couple/cfd in LIGGGHTS!
+
+Some warnings may be given for particles that have not been located - this is due to LIGGGHTS' treatment of domain crossers. 
+
+[Related commands:]
+
+"dataExchangeModel"_dataExchangeModel.html
+

doc/forceModel_gradPForceSmooth.txt

@@ -0,0 +1,77 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+forceModel gradPForceSmooth command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+forceModels
+(
+    gradPForceSmooth;
+);
+gradPForceSmoothProps
+\{
+    pFieldName         "pressure";
+    velocityFieldName  "U";
+    useAddedMass       scalar1;
+    treatForceExplicit switch1;
+    treatForceDEM      switch2;
+    interpolation      switch3;
+    smoothingModel     "smoothingModel";
+\} :pre
+
+{pressure} = name of the finite volume fluid pressure field :ulb,l
+{U} = name of the finite volume fluid velocity field :l
+{scalar1} = (optional, default 0) coefficient of added mass accounted for :l
+{switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{switch3} = (optional, default false) flag to use interpolated pressure values :l
+{smoothingModel} = name of smoothing model :l
+:ule
+
+[Examples:]
+
+forceModels
+(
+    gradPForceSmooth;
+);
+gradPForceSmoothProps
+\{
+    pFieldName        "p_rgh";
+    velocityFieldName "U";
+    interpolation     false;
+    smoothingModel    "temporalSmoothing";
+    temporalSmoothingProps
+    \{
+        lowerLimit 0.1;
+        upperLimit 1e10;
+        refField   "p_rgh";
+        gamma      1.0;
+    \} :pre
+
+\} :pre
+[Description:]
+
+The {gradPForceSmooth} model calculates the particle based pressure gradient
+force identically to the "gradPForce"_forceModel_gradPForce.html model but
+allows smoothing of the pressure prior to the force calculation (without
+altering the original pressure field).
+Any smoothing model can be used and does not have to be the same as specified in
+couplingProperties. Properties for the smoothing model have to be specified in a
+sub-dictionary within {gradPForceSmoothProps}.
+
+[Restrictions:]
+
+A volScalarField "pSmooth" MUST be specified in the initial time directory!
+
+[Related commands:]
+
+"forceModel"_forceModel.html, "gradPForce"_forceModel_gradPForce.html
+

doc/forceModel_surfaceTensionForce.txt

@@ -0,0 +1,55 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+forceModel surfaceTensionForce command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+forceModels
+(
+    surfaceTensionForce;
+);
+surfaceTensionForceProps
+\{
+    stfFieldName "surfaceTensionField";
+\} :pre
+
+{surfaceTensionField} = name of the surface tension force field :ulb,l
+:ule
+
+[Examples:]
+
+forceModels
+(
+    surfaceTensionForce;
+);
+surfaceTensionForceProps
+\{
+    stfFieldName "surfaceTensionForce";
+\} :pre
+
+[Description:]
+
+The force model calculates the surface tension force acting on each DEM particle
+based on a pre-calculated surface tension force as V_particle * F^sigma. When
+used in conjunction with the "cfdemSolverMultiphase"_cfdemSolverMultiphase.html
+solver, the surface tension force is calculated with the CSF (continuum surface
+force) model (see Brackbill et al. (1992): "A continuum method for modeling
+surface tension", J. Comput. Phys.).
+
+[Restrictions:]
+
+Has to be used with a multiphase solver that calculates the surface tension
+force, e.g. {cfdemSolverMultiphase}.
+
+[Related commands:]
+
+"forceModel"_forceModel.html, "cfdemSolverMultiphase"_cfdemSolverMultiphase.html
+

doc/smoothingModel_temporalSmoothing.txt

@@ -0,0 +1,62 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+smoothingModel temporalSmoothing command :h3
+
+[Syntax:]
+
+Defined in dictionary depending on the application.
+
+smoothingModel temporalSmoothing;
+temporalSmoothingProps
+\{
+    lowerLimit number1;
+    upperLimit number2;
+    refField   referenceField;
+    gamma      smoothingStrength;
+\} :pre
+
+{number1} = scalar fields will be bound to this lower value :ulb,l
+{number2} = scalar fields will be bound to this upper value :l
+{referenceField} = reference to the un-smoothed field required for the relaxation operation :l
+{smoothingStrength} = control parameter for the smoothing, lower value yields stronger smoothing (gamma = 1 results in an equal contribution from the un-smoothed and smoothed fields)  :l
+:ule
+
+[Examples:]
+
+temporalSmoothingProps
+\{
+    lowerLimit     0.1;
+    upperLimit     1e10;
+    referenceField "p";
+    gamma          1.0;
+\} :pre
+
+[Description:]
+
+The {temporalSmoothing} model is a smoothing model that utilizes temporal
+relaxation of a desired quantity. This model can be used to filter out
+high-frequency fluctuations (e.g. numerical noise) controlled via the control
+parameter gamma.
+Note that this model does NOT smooth the calculated fields, instead smoothing is
+performed on a separate (smooth) field which uses the calculated (un-smooth)
+field as a reference.
+Thus its usage is limited and CANNOT be used to smooth the exchange fields
+similar to other smoothing models.
+For further information see Vångö et al., "Unresolved CFD-DEM modeling of
+multiphase flow in densely packed particle beds", Appl. Math. Model. (2018).
+
+[Restrictions:]
+
+This model does NOT smooth the calculated fields and can therefore NOT be used
+as a general smoothing model to smoothen the exchange fields.
+Attempting this will generate an error.
+
+[Related commands:]
+
+"smoothingModel"_smoothingModel.html
+

etc/bashrc

@@ -17,7 +17,7 @@
 #------------------------------------------------------------------------------
 
 export CFDEM_PROJECT=CFDEM
-export CFDEM_VERSION=18.10
+export CFDEM_VERSION=19.02
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade

etc/cshrc

@@ -15,7 +15,7 @@
 #------------------------------------------------------------------------------
 
 setenv CFDEM_PROJECT CFDEM
-setenv CFDEM_VERSION 18.10
+setenv CFDEM_VERSION 19.02
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade

etc/library-list.txt

@@ -1,3 +1,4 @@
 lagrangian/cfdemParticle/dir
 lagrangian/cfdemParticleComp/dir
 finiteVolume/dir
+../applications/solvers/cfdemSolverMultiphase/multiphaseMixture/dir

etc/solver-list.txt

@@ -5,3 +5,4 @@ cfdemSolverRhoSimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
 cfdemSolverRhoPimpleChem/dir
+cfdemSolverMultiphase/dir

src/lagrangian/cfdemParticle/Make/files

@@ -80,6 +80,8 @@ $(forceModels)/Fines/FanningDynFines.C
 $(forceModels)/Fines/ErgunStatFines.C
 $(forceModels)/granKineticEnergy/granKineticEnergy.C
 $(forceModels)/pdCorrelation/pdCorrelation.C
+$(forceModels)/surfaceTensionForce/surfaceTensionForce.C
+$(forceModels)/gradPForceSmooth/gradPForceSmooth.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
@@ -147,6 +149,8 @@ $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
 $(dataExchangeModels)/twoWayFiles/twoWayFiles.C
 $(dataExchangeModels)/noDataExchange/noDataExchange.C
 $(dataExchangeModels)/twoWayMPI/twoWayMPI.C
+$(dataExchangeModels)/twoWayOne2One/twoWayOne2One.C
+$(dataExchangeModels)/twoWayOne2One/one2one.C
 
 $(averagingModels)/averagingModel/averagingModel.C
 $(averagingModels)/averagingModel/newAveragingModel.C
@@ -169,5 +173,6 @@ $(smoothingModels)/smoothingModel/smoothingModel.C
 $(smoothingModels)/smoothingModel/newSmoothingModel.C
 $(smoothingModels)/noSmoothing/noSmoothing.C
 $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C
+$(smoothingModels)/temporalSmoothing/temporalSmoothing.C
 
 LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME)

src/lagrangian/cfdemParticle/Make/options

@@ -34,6 +34,4 @@ LIB_LIBS = \
     -lmpi_cxx \
     -Wl,-rpath,$(CFDEM_LIGGGHTS_BIN_DIR) \
     -L$(CFDEM_LIGGGHTS_BIN_DIR) \
-    -lliggghts \
-    -L$(CFDEM_Many2ManyLIB_PATH) \
-    -lcoupleMany2Many
+    -lliggghts

src/lagrangian/cfdemParticle/cfdTools/checkModelType.H

@@ -21,7 +21,7 @@
         found=false;
         forAll(particleCloud.forceModels(),i)
         {
-            if(particleCloud.forceModels()[i]=="gradPForce")
+            if(particleCloud.forceModels()[i]=="gradPForce" || particleCloud.forceModels()[i]=="gradPForceSmooth")
                 found=true;
         }
         if(!found)
@@ -56,7 +56,7 @@
         found=false;
         forAll(particleCloud.forceModels(),i)
         {
-            if(particleCloud.forceModels()[i]=="gradPForce" || particleCloud.forceModels()[i]=="viscForce")
+            if(particleCloud.forceModels()[i]=="gradPForce" || particleCloud.forceModels()[i]=="gradPForceSmooth" || particleCloud.forceModels()[i]=="viscForce")
                 found=true;
         }
         if(found)
@@ -80,7 +80,7 @@
         found=false;
         forAll(particleCloud.forceModels(),i)
         {
-            if(particleCloud.forceModels()[i]=="gradPForce")
+            if(particleCloud.forceModels()[i]=="gradPForce" || particleCloud.forceModels()[i]=="gradPForceSmooth")
                 found=true;
         }
         if(!found)
@@ -99,3 +99,6 @@
         Warning << "You chose model type -none- you might get erroneous results!" << endl;
     else
         FatalError <<"no suitable model type specified:" << modelType << "\n" << abort(FatalError);
+
+    if (particleCloud.smoothingM().type() == "temporalSmoothing")
+        FatalError << "the temporalSmoothing model does not support smoothing of the exchange fields, please see documentation!" << endl;

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -34,8 +34,8 @@ Description
 #ifndef versionInfo_H
 #define versionInfo_H
 
-word CFDEMversion="PFM 18.10";
-word compatibleLIGGGHTSversion="PFM 18.10";
+word CFDEMversion="PFM 19.02";
+word compatibleLIGGGHTSversion="PFM 19.02";
 word OFversion="4.x";
 
 Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;

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

@@ -312,7 +312,7 @@ void diffusionCoefficient::execute()
                                 TotalFraction_[i]   +=  Xfluid_[j]/dBinary_;
 
                                 // dCoeff -- diffusion component of diffusant gas
-                                MixtureBinaryDiffusion_[i]  =   (1.0-XfluidDiffusant_[i])/TotalFraction_[i];                             
+                                MixtureBinaryDiffusion_[i]  =   (1.0-XfluidDiffusant_[i])/TotalFraction_[i];
 
                                 if(verbose_)
                                 {

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

@@ -72,7 +72,7 @@ private:
     // gas pressure at particle location
     word pressureFieldName_;
 
-    const volScalarField& P_; 
+    const volScalarField& P_;
 
     word partPressureName_;
 

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

@@ -109,8 +109,6 @@ void reactantPerParticle::reAllocMyArrays() const
 
 void reactantPerParticle::execute()
 {
-   
-    
     loopCounter_++;
     if (loopCounter_ % Nevery_ != 0)
     {
@@ -121,7 +119,7 @@ void reactantPerParticle::execute()
 
     particlesPerCell_ *= 0.0;
 
-    label  cellI=0; 
+    label  cellI=0;
     scalar voidfraction(1);
     scalar cellvolume(0.0);
     scalar particlesPerCell(1.0);
@@ -133,7 +131,7 @@ void reactantPerParticle::execute()
         if (cellI >= 0)
         {
             particlesPerCell_[cellI] += 1.0;
-	}
+        }
     }
 
     // no fill array and communicate it
@@ -151,7 +149,7 @@ void reactantPerParticle::execute()
 
         // give DEM data
         particleCloud_.dataExchangeM().giveData("reactantPerParticle", "scalar-atom", reactantPerParticle_);
-       
+
         Info << "give data done" << endl;
 }
 

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

@@ -245,17 +245,25 @@ public:
         virtual int getNumberOfTypes() const;
         virtual double* getTypeVol() const;
 
-        inline void setPositions(label n,double* pos) const
+        inline void setPositions(label n,double* pos)
         {
             for (int i=0;i<n;i++)
                 for (int j=0;j<3;j++)
                     particleCloud_.positions_[i][j]=pos[i*3+j];
         }
-        inline void setCellIDs(label n,int* ID) const
+
+        inline void setCellIDs(label n,int* ID)
         {
             for (int i=0;i<n;i++)
                 particleCloud_.cellIDs_[i][0]=ID[i];
         }
+        inline void setCellIDs(labelList const& ids)
+        {
+            for (int i = 0; i < ids.size(); i++)
+            {
+                particleCloud_.cellIDs_[i][0] = ids[i];
+            }
+        }
 
 
         virtual scalar getCG() const { Warning << "getCG() not executed correctly!" << endl; return 1.; }

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

@@ -0,0 +1,286 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+Copyright 2018-     Paul Kieckhefen, TUHH
+
+\*---------------------------------------------------------------------------*/
+
+//#define O2ODEBUG
+
+#ifdef O2ODEBUG
+#include <iostream>
+#endif
+
+#include <mpi.h>
+#include "one2one.H"
+
+template <typename T>
+struct mpi_type_wrapper {
+  MPI_Datatype mpi_type;
+  mpi_type_wrapper();
+};
+template <> mpi_type_wrapper<float>::mpi_type_wrapper()
+: mpi_type(MPI_FLOAT) {}
+template <> mpi_type_wrapper<double>::mpi_type_wrapper()
+: mpi_type(MPI_DOUBLE) {}
+template <> mpi_type_wrapper<int>::mpi_type_wrapper()
+: mpi_type(MPI_INT) {}
+
+/* ---------------------------------------------------------------------- */
+
+One2One::One2One(MPI_Comm caller)
+:
+ncollected_(-1),
+comm_(caller),
+nsrc_procs_(-1),
+src_procs_(nullptr),
+ndst_procs_(-1),
+dst_procs_(nullptr),
+nlocal_(-1),
+natoms_(nullptr),
+request_(nullptr),
+status_(nullptr)
+{
+  MPI_Comm_rank(comm_,&me_);
+  MPI_Comm_size(comm_,&nprocs_);
+}
+
+/* ---------------------------------------------------------------------- */
+
+One2One::~One2One()
+{
+  deallocate();
+}
+
+/* ----------------------------------------------------------------------
+  communicate particle ids based on processor communication pattern
+------------------------------------------------------------------------- */
+
+void One2One::setup
+(
+  int nsrc_procs,
+  int *src_procs,
+  int ndst_procs,
+  int *dst_procs,
+  int nlocal
+)
+{
+  // free any previous one2one info
+  deallocate();
+
+  src_procs_ = src_procs;
+  nsrc_procs_ = nsrc_procs;
+  dst_procs_ = dst_procs;
+  ndst_procs_ = ndst_procs;
+  nlocal_ = nlocal;
+
+  // gather number of ids for reserving memory
+  natoms_ = new int[nprocs_];
+  MPI_Allgather // may be replaced by send/irecv
+  (
+    &nlocal_,
+    1,
+    MPI_INT,
+    natoms_,
+    1,
+    MPI_INT,
+    comm_
+  );
+
+  ncollected_ = 0;
+  int nrequests = 0;
+  for (int i = 0; i < nsrc_procs_; i++)
+  {
+    if (natoms_[src_procs_[i]] > 0)
+    {
+      ncollected_ += natoms_[src_procs_[i]];
+
+      if (src_procs_[i] != me_) // no receive for on-proc info
+      {
+        nrequests++;
+      }
+    }
+  }
+
+  if (nrequests > 0)
+  {
+    request_ = new MPI_Request[nrequests];
+    status_ = new MPI_Status[nrequests];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   src: what is present on this proc
+   dst: what is received from other procs
+   all comm according to map set up in setup(...)
+------------------------------------------------------------------------- */
+template  <typename T>
+void One2One::exchange(T *&src, T *&dst, int data_length)
+{
+  mpi_type_wrapper<T> wrap;
+
+  // post receives
+  int offset_local = -1;
+  int offset = 0;
+  int requesti = 0;
+  for (int i = 0; i < nsrc_procs_; i++)
+  {
+    // do post a receives for procs who own particles
+    if (natoms_[src_procs_[i]] > 0)
+    {
+      if (src_procs_[i] != me_)
+      {
+  #ifdef O2ODEBUG
+  std::cout<< "[" << me_ << "]"
+           << " RCV " << i
+           << " of "  << nsrc_procs_
+           << " from: " << src_procs_[i]
+           << " natoms_[src_procs_[i]] " << natoms_[src_procs_[i]]
+           << " datalength " << data_length
+           << " offset " << offset
+           << std::endl;
+  #endif
+        MPI_Irecv
+        (
+          &dst[offset],
+          natoms_[src_procs_[i]]*data_length,
+          wrap.mpi_type,
+          src_procs_[i],
+          MPI_ANY_TAG,
+          comm_,
+          &request_[requesti]
+        );
+        requesti++;
+      }
+      else // data is available on-proc
+      {
+        offset_local = offset;
+      }
+    }
+    offset += natoms_[src_procs_[i]]*data_length;
+  }
+
+  // make sure all receives are posted
+  MPI_Barrier(comm_);
+
+  // blocking sends - do nonblocking instead
+  //                  since doing many-2-many here?
+  // only do sends if I have particles
+  if (nlocal_ > 0)
+  {
+    for (int i = 0; i < ndst_procs_; i++)
+    {
+      if (dst_procs_[i] != me_)
+      {
+    #ifdef O2ODEBUG
+    std::cout<< "[" << me_ << "]"
+             << " SEND to: " << dst_procs_[i]
+             << " nlocal_ " << nlocal_
+             << " data_length " << data_length
+             << std::endl;
+    #endif
+        MPI_Send
+        (
+          src,
+          nlocal_*data_length,
+          wrap.mpi_type,
+          dst_procs_[i],
+          0,
+          comm_
+        );
+      }
+    }
+  }
+
+  // only wait if requests were actually posted
+  if (requesti > 0)
+    MPI_Waitall(requesti, request_, status_);
+
+  // copy on-proc data
+  if (offset_local > -1)
+  {
+    const int max_locali = nlocal_ * data_length;
+    for
+    (
+      int locali = 0;
+      locali < max_locali;
+      locali++
+    )
+    {
+      dst[locali+offset_local] = src[locali];
+    }
+  }
+}
+
+template void One2One::exchange<int>(int*&, int*&, int);
+template void One2One::exchange<double>(double*&, double*&, int);
+
+// there should be a way to do this without copying data
+template  <typename T>
+void One2One::exchange(T **&src, T **&dst, int data_length)
+{
+  mpi_type_wrapper<T> wrap;
+
+  T* tmp_dst = new T[ncollected_*data_length];
+  T* tmp_src = new T[nlocal_*data_length];
+
+  for (int i = 0; i < nlocal_; i++)
+    for (int j = 0; j < data_length; j++)
+      tmp_src[data_length*i+j] = src[i][j];
+
+  exchange<T>(tmp_src, tmp_dst, data_length);
+
+  for (int i = 0; i < ncollected_; i++)
+    for (int j = 0; j < data_length; j++)
+      dst[i][j] = tmp_dst[data_length*i+j];
+
+  delete [] tmp_src;
+  delete [] tmp_dst;
+}
+template void One2One::exchange<int>(int**&, int**&, int);
+template void One2One::exchange<double>(double**&, double**&, int);
+
+
+template  <typename T>
+void One2One::exchange(T **&src, T *&dst, int data_length)
+{
+  mpi_type_wrapper<T> wrap;
+
+  T* tmp_src = new T[nlocal_*data_length];
+
+  for (int i = 0; i < nlocal_; i++)
+    for (int j = 0; j < data_length; j++)
+      tmp_src[data_length*i+j] = src[i][j];
+
+  exchange<T>(tmp_src, dst, data_length);
+
+
+  delete [] tmp_src;
+}
+template void One2One::exchange<int>(int**&, int*&, int);
+template void One2One::exchange<double>(double**&, double*&, int);
+
+/* ---------------------------------------------------------------------- */
+
+void One2One::deallocate()
+{
+  delete [] src_procs_;
+  delete [] dst_procs_;
+  delete [] natoms_;
+
+  delete [] request_;
+  delete [] status_;
+}

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

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+Copyright 2018-     Paul Kieckhefen, TUHH
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ONE2ONE_H
+#define ONE2ONE_H
+
+#include <mpi.h>
+
+class One2One
+{
+ public:
+  One2One(MPI_Comm);
+  ~One2One();
+
+  void setup
+  (
+    int nsrc_procs,
+    int *src_procs,
+    int ndst_procs,
+    int* dst_procs,
+    int nlocal
+  );
+
+  template <typename T>
+  void exchange(T *&, T *&, int data_length=1);
+  template <typename T>
+  void exchange(T **&, T **&, int data_length=1);
+  template <typename T>
+  void exchange(T **&, T *&, int data_length=1);
+
+  int ncollected_;            // # of ids in from group
+
+ protected:
+
+  MPI_Comm comm_;
+  int me_, nprocs_;                // rank and size
+
+  // communication partners
+  int  nsrc_procs_;                // # of off-processor IDs
+  int* src_procs_;                 // procs I receive data from
+  int  ndst_procs_;                // # of off-processor IDs
+  int* dst_procs_;                 // procs I receive data from
+
+  int nlocal_;                     // # particle ids I own
+  int* natoms_;
+
+  MPI_Request* request_;
+  MPI_Status* status_;
+
+  void deallocate();
+};
+
+#endif

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

@@ -0,0 +1,941 @@
+/*---------------------------------------------------------------------------*\
+    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(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+Contributing authors
+    Paul Kieckhefen (TUHH)  2018-
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "twoWayOne2One.H"
+#include "addToRunTimeSelectionTable.H"
+#include "clockModel.H"
+#include "pair.h"
+#include "force.h"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(twoWayOne2One, 0);
+
+addToRunTimeSelectionTable
+(
+    dataExchangeModel,
+    twoWayOne2One,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+twoWayOne2One::twoWayOne2One
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    dataExchangeModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    thisLigPartner_(0),
+    thisFoamPartner_(0),
+    lig2foam_(nullptr),
+    foam2lig_(nullptr),
+    lig2foam_mask_(nullptr),
+    lig2foam_ids_(nullptr),
+    foam2lig_ids_(nullptr),
+    lig2foam_vec_tmp_(nullptr),
+    lig2foam_scl_tmp_(nullptr),
+    foam2lig_vec_tmp_(nullptr),
+    foam2lig_scl_tmp_(nullptr),
+    staticProcMap_(propsDict_.lookupOrDefault<Switch>("useStaticProcMap", false)),
+    cellIdComm_(propsDict_.lookupOrDefault<Switch>("useCellIdComm", false)),
+    my_prev_cell_ids_fix_(nullptr),
+    verbose_(propsDict_.lookupOrDefault("verbose", false)),
+    lmp(nullptr)
+{
+    Info<<"Starting up LIGGGHTS for first time execution"<<endl;
+
+    comm_liggghts_ = MPI_COMM_WORLD;
+
+    // read path from dictionary
+    const fileName liggghtsPath(propsDict_.lookup("liggghtsPath"));
+
+    // open LIGGGHTS input script
+    Info<<"Executing input script '"<< liggghtsPath.c_str() <<"'"<<endl;
+    lmp = new LAMMPS_NS::LAMMPS(0,nullptr,comm_liggghts_);
+    lmp->input->file(liggghtsPath.c_str());
+
+    // get DEM time step size
+    DEMts_ = lmp->update->dt;
+    checkTSsize();
+
+    // calculate boundingBox of FOAM subdomain
+    primitivePatch tmpBoundaryFaces
+    (
+        SubList<face>
+        (
+            sm.mesh().faces(),
+            sm.mesh().nFaces() - sm.mesh().nInternalFaces(),
+            sm.mesh().nInternalFaces()
+        ),
+        sm.mesh().points()
+    );
+    typedef PrimitivePatch<face, List, const pointField> bPatch;
+    bPatch boundaryFaces
+    (
+        tmpBoundaryFaces.localFaces(),
+        tmpBoundaryFaces.localPoints()
+    );
+    thisFoamBox_ = treeBoundBox(boundaryFaces.localPoints());
+    if (staticProcMap_)
+    {
+        createProcMap();
+    }
+
+    if (cellIdComm_)
+    {
+        my_prev_cell_ids_fix_ = static_cast<LAMMPS_NS::FixPropertyAtom*>
+        (   lmp->modify->find_fix_property
+            (
+                "prev_cell_ids",
+                "property/atom",
+                "scalar",
+                0,
+                0,
+                "cfd coupling",
+                true
+            )
+        );
+    }
+}
+
+void twoWayOne2One::createProcMap()
+{
+    List<treeBoundBox> foamBoxes(Pstream::nProcs());
+    foamBoxes[Pstream::myProcNo()] = thisFoamBox_;
+    Pstream::gatherList(foamBoxes);
+    Pstream::scatterList(foamBoxes);
+
+    // calculate bounding box of LIG subdomain
+    // this may have to move to couple when dynamic LB occurs
+    List<boundBox> ligBoxes(Pstream::nProcs());
+    double** ligbb = o2o_liggghts_get_boundingbox(lmp);
+    boundBox thisLigBox
+    (
+        point(ligbb[0][0], ligbb[0][1], ligbb[0][2]),
+        point(ligbb[1][0], ligbb[1][1], ligbb[1][2])
+    );
+    ligBoxes[Pstream::myProcNo()] = thisLigBox;
+    Pstream::gatherList(ligBoxes);
+    Pstream::scatterList(ligBoxes);
+
+    thisLigPartner_.clear();
+    thisFoamPartner_.clear();
+
+    // detect LIG subdomains which this FOAM has to interact with
+    forAll(ligBoxes, ligproci)
+    {
+        if (thisFoamBox_.overlaps(ligBoxes[ligproci]))
+        {
+            thisLigPartner_.append(ligproci);
+        }
+    }
+    // detect FOAM subdomains this LIG has to interact with
+    // TODO: refactor to invert this list here
+    forAll(foamBoxes, foamproci)
+    {
+        if (thisLigBox.overlaps(foamBoxes[foamproci]))
+        {
+            thisFoamPartner_.append(foamproci);
+        }
+    }
+
+    if (verbose_)
+    {
+        Pout<< "FOAM bounding box: " << thisFoamBox_
+            << " LIG bounding box: " << thisLigBox
+            << nl
+            << "FOAM comm partners: " << thisFoamPartner_
+            << " LIG comm partners: " << thisLigPartner_
+            << endl;
+    }
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+twoWayOne2One::~twoWayOne2One()
+{
+    delete foam2lig_;
+    delete lig2foam_;
+
+    destroy(lig2foam_ids_);
+    destroy(foam2lig_ids_);
+
+    destroy(lig2foam_vec_tmp_);
+    destroy(lig2foam_scl_tmp_);
+    destroy(foam2lig_vec_tmp_);
+    destroy(foam2lig_scl_tmp_);
+
+    delete lmp;
+}
+
+// * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
+void twoWayOne2One::getData
+(
+    word name,
+    word type,
+    double ** const& field,
+    label /*step*/
+) const
+{
+    if (name == "x") // the location is transferred by couple()
+    {
+        return;
+    }
+    if (type == "vector-atom")
+    {
+        double **tmp_= static_cast<double **>(lammps_extract_atom(lmp,name.c_str()));
+        if (!tmp_)
+        {
+            LAMMPS_NS::Fix *fix = nullptr;
+            fix = lmp->modify->find_fix_property
+            (
+                name.c_str(),
+                "property/atom",
+                "vector",
+                0,
+                0,
+                "cfd coupling",
+                false
+            );
+            if (fix)
+            {
+                tmp_ = static_cast<LAMMPS_NS::FixPropertyAtom*>(fix)->array_atom;
+            }
+            else
+            {
+                Warning<< "coupling fix not found!" << endl;
+            }
+
+            if (!tmp_)
+            {
+                FatalError<< "find_fix_property " << name
+                          << " array_atom not found."
+                          << abort(FatalError);
+            }
+        }
+
+        lig2foam_->exchange<double>
+        (
+            tmp_,
+            lig2foam_vec_tmp_,
+            3
+        );
+        extractCollected<double>
+        (
+            lig2foam_vec_tmp_,
+            const_cast<double**&>(field),
+            3
+        );
+    }
+    else if (type == "scalar-atom")
+    {
+        double *tmp_ = static_cast<double *>(lammps_extract_atom(lmp,name.c_str()));
+        if (!tmp_)
+        {
+            LAMMPS_NS::Fix *fix = nullptr;
+            fix = lmp->modify->find_fix_property
+            (
+                name.c_str(),
+                "property/atom",
+                "scalar",
+                0,
+                0,
+                "cfd coupling",
+                true
+            );
+
+            if (fix)
+            {
+                tmp_ = static_cast<LAMMPS_NS::FixPropertyAtom*>(fix)->vector_atom;
+            }
+            else
+            {
+                FatalError<< "coupling fix not found!" << abort(FatalError);
+            }
+
+            if (!tmp_)
+            {
+                FatalError<< "find_fix_property " << name
+                          << " vector_atom not found."
+                          << abort(FatalError);
+            }
+        }
+        lig2foam_->exchange<double>
+        (
+            tmp_,
+            lig2foam_scl_tmp_
+        );
+        extractCollected<double>
+        (
+            lig2foam_scl_tmp_,
+            const_cast<double**&>(field)
+        );
+    }
+    else
+    {
+        FatalError << "requesting type " << type << " and name " << name << abort(FatalError);
+    }
+}
+
+void twoWayOne2One::getData
+(
+    word name,
+    word type,
+    int ** const& field,
+    label /*step*/
+) const
+{
+    FatalError << "do not use this getData!!!" << abort(FatalError);
+/*
+    o2o_data_liggghts_to_of
+    (
+        name.c_str(),
+        type.c_str(),
+        lmp,
+        (void*&) field,
+        "int",
+        comm_liggghts_
+    );
+*/
+}
+
+void twoWayOne2One::giveData
+(
+    word name,
+    word type,
+    double ** const& field,
+    const char* datatype
+) const
+{
+    if (type == "vector-atom")
+    {
+        foam2lig_->exchange
+        (
+            const_cast<double**&>(field),
+            foam2lig_vec_tmp_,
+            3
+        );
+        o2o_data_of_to_liggghts
+        (
+            name.c_str(),
+            type.c_str(),
+            lmp,
+            foam2lig_vec_tmp_,
+            datatype,
+            foam2lig_ids_,
+            foam2lig_->ncollected_
+        );
+    }
+    else if (type == "scalar-atom")
+    {
+        foam2lig_->exchange
+        (
+            const_cast<double**&>(field),
+            foam2lig_scl_tmp_,
+            1
+        );
+        o2o_data_of_to_liggghts
+        (
+            name.c_str(),
+            type.c_str(),
+            lmp,
+            foam2lig_scl_tmp_,
+            datatype,
+            foam2lig_ids_,
+            foam2lig_->ncollected_
+        );
+    }
+    else
+    {
+        FatalError<< "twoWayMany2Many::giveData requested type " << type
+                  << " not implemented!"
+                  << abort(FatalError);
+    }
+}
+
+//============
+// double **
+void twoWayOne2One::allocateArray
+(
+    double**& array,
+    double initVal,
+    int width,
+    int length
+) const
+{
+    int len = max(length,1);
+    lmp->memory->grow(array, len, width, "o2o:dbl**");
+    for (int i = 0; i < len; i++)
+        for (int j = 0; j < width; j++)
+            array[i][j] = initVal;
+}
+
+void twoWayOne2One::allocateArray
+(
+    double**& array,
+    double initVal,
+    int width,
+    const char* length
+) const
+{
+    int len = max(particleCloud_.numberOfParticles(),1);
+    lmp->memory->grow(array, len, width, "o2o:dbl**:autolen");
+    for (int i = 0; i < len; i++)
+        for (int j = 0; j < width; j++)
+            array[i][j] = initVal;
+}
+
+void inline twoWayOne2One::destroy(double** array,int len) const
+{
+    lmp->memory->destroy(array);
+}
+
+//============
+// int **
+void twoWayOne2One::allocateArray
+(
+    int**& array,
+    int initVal,
+    int width,
+    int length
+) const
+{
+    int len = max(length,1);
+    lmp->memory->grow(array, len, width, "o2o:int**");
+    for (int i = 0; i < len; i++)
+        for (int j = 0; j < width; j++)
+            array[i][j] = initVal;
+}
+
+void twoWayOne2One::allocateArray
+(
+    int**& array,
+    int initVal,
+    int width,
+    const char* length
+) const
+{
+    int len = max(particleCloud_.numberOfParticles(),1);
+    lmp->memory->grow(array, len, width, "o2o:int**:autolen");
+    for (int i = 0; i < len; i++)
+        for (int j = 0; j < width; j++)
+            array[i][j] = initVal;
+}
+
+void inline twoWayOne2One::destroy(int** array,int len) const
+{
+    lmp->memory->destroy(array);
+}
+
+//============
+// double *
+void twoWayOne2One::allocateArray(double*& array, double initVal, int length) const
+{
+    int len = max(length,1);
+    lmp->memory->grow(array, len, "o2o:dbl*");
+    for (int i = 0; i < len; i++)
+        array[i] = initVal;
+}
+
+void inline twoWayOne2One::destroy(double* array) const
+{
+    lmp->memory->destroy(array);
+}
+
+//==============
+// int *
+void twoWayOne2One::allocateArray(int*& array, int initVal, int length) const
+{
+    int len = max(length,1);
+    lmp->memory->grow(array, len, "o2o:int*");
+    for (int i = 0; i < len; i++)
+        array[i] = initVal;
+}
+
+void inline twoWayOne2One::destroy(int* array) const
+{
+    lmp->memory->destroy(array);
+}
+//==============
+
+bool twoWayOne2One::couple(int i)
+{
+    bool coupleNow = false;
+    if (i==0)
+    {
+        couplingStep_++;
+        coupleNow = true;
+
+
+        // run commands from liggghtsCommands dict
+        Info<< "Starting up LIGGGHTS" << endl;
+        particleCloud_.clockM().start(3,"LIGGGHTS");
+
+        // check if liggghtsCommandModels with exaxt timing are being run
+        bool exactTiming(false);
+        int runComNr = -10;
+        DynamicList<scalar> interruptTimes(0);
+        DynamicList<int> DEMstepsToInterrupt(0);
+        DynamicList<int> lcModel(0);
+
+        forAll(particleCloud_.liggghtsCommandModelList(),i)
+        {
+            // Check if exact timing is needed
+            // get time for execution
+            // store time for execution in list
+            if(particleCloud_.liggghtsCommand(i).exactTiming())
+            {
+                exactTiming = true;
+                DynamicList<scalar> h
+                  = particleCloud_.liggghtsCommand(i).executionsWithinPeriod
+                (
+                    TSstart(),
+                    TSend()
+                );
+
+                forAll(h,j)
+                {
+                    // save interrupt times (is this necessary)
+                    interruptTimes.append(h[j]);
+
+                    // calc stepsToInterrupt
+                    DEMstepsToInterrupt.append(DEMstepsTillT(h[j]));
+
+                    // remember which liggghtsCommandModel to run
+                    lcModel.append(i);
+                }
+
+                // make cumulative
+                label len = DEMstepsToInterrupt.size();
+                label ind(0);
+                forAll(DEMstepsToInterrupt,i)
+                {
+                    ind = len - i - 1;
+                    if(ind > 0)
+                    {
+                        DEMstepsToInterrupt[ind] -= DEMstepsToInterrupt[ind-1];
+                    }
+                }
+
+                Info<< "Foam::twoWayOne2One::couple(i): interruptTimes="      << interruptTimes      << nl
+                    << "Foam::twoWayOne2One::couple(i): DEMstepsToInterrupt=" << DEMstepsToInterrupt << nl
+                    << "Foam::twoWayOne2One::couple(i): lcModel="             << lcModel
+                    << endl;
+            }
+
+            if(particleCloud_.liggghtsCommand(i).type() == "runLiggghts")
+            {
+                runComNr = i;
+            }
+        }
+
+        // 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(0);
+                Info<< "Executing run command: '"<< command <<"'"<< endl;
+                lmp->input->one(command);
+
+                // run liggghts command with exact timing
+                command = particleCloud_.liggghtsCommand(lcModel[j]).command(0);
+                Info << "Executing command: '"<< command <<"'"<< endl;
+                lmp->input->one(command);
+            }
+
+            // do the run
+            if(particleCloud_.liggghtsCommand(runComNr).runCommand(couplingStep()))
+            {
+                particleCloud_.liggghtsCommand(runComNr).set(couplingInterval() - DEMstepsRun);
+                const char* command = particleCloud_.liggghtsCommand(runComNr).command(0);
+                Info<< "Executing run command: '"<< command <<"'"<< endl;
+                lmp->input->one(command);
+            }
+
+            // do the other non exact timing models
+            forAll(particleCloud_.liggghtsCommandModelList(),i)
+            {
+                if
+                (
+                  ! particleCloud_.liggghtsCommand(i).exactTiming() &&
+                    particleCloud_.liggghtsCommand(i).runCommand(couplingStep())
+                )
+                {
+                    commandLines=particleCloud_.liggghtsCommand(i).commandLines();
+                    for(int j=0;j<commandLines;j++)
+                    {
+                        const char* command = particleCloud_.liggghtsCommand(i).command(j);
+                        Info << "Executing command: '"<< command <<"'"<< endl;
+                        lmp->input->one(command);
+                    }
+                }
+            }
+        }
+        // no model with exact timing exists
+        else
+        {
+            forAll(particleCloud_.liggghtsCommandModelList(),i)
+            {
+                if(particleCloud_.liggghtsCommand(i).runCommand(couplingStep()))
+                {
+                    commandLines=particleCloud_.liggghtsCommand(i).commandLines();
+                    for(int j=0;j<commandLines;j++)
+                    {
+                        const char* command = particleCloud_.liggghtsCommand(i).command(j);
+                        Info << "Executing command: '"<< command <<"'"<< endl;
+                        lmp->input->one(command);
+                    }
+                }
+            }
+        }
+
+        particleCloud_.clockM().stop("LIGGGHTS");
+        Info<< "LIGGGHTS finished" << endl;
+
+        if (!staticProcMap_)
+        {
+            createProcMap();
+        }
+
+        setupLig2FoamCommunication();
+
+        locateParticles();
+
+        setupFoam2LigCommunication();
+
+        if (verbose_)
+        {
+            Pout<< "FOAM owns " << getNumberOfParticles()
+                << " LIG owns " << lmp->atom->nlocal
+                << nl
+                << "FOAM collects " << lig2foam_->ncollected_
+                << " LIG collects " << foam2lig_->ncollected_
+                << endl;
+        }
+    }
+
+    return coupleNow;
+}
+
+void twoWayOne2One::setupLig2FoamCommunication()
+{
+    int* src_procs = new int[thisLigPartner_.size()];
+    for (int proci = 0; proci < thisLigPartner_.size(); proci++)
+    {
+        src_procs[proci] = thisLigPartner_[proci];
+    }
+    int* dst_procs = new int[thisFoamPartner_.size()];
+    for (int proci = 0; proci < thisFoamPartner_.size(); proci++)
+    {
+        dst_procs[proci] = thisFoamPartner_[proci];
+    }
+
+    delete lig2foam_;
+    lig2foam_ = new One2One(comm_liggghts_);
+    lig2foam_->setup
+    (
+      thisLigPartner_.size(),
+      src_procs,
+      thisFoamPartner_.size(),
+      dst_procs,
+      lmp->atom->nlocal
+    );
+    allocateArray
+    (
+        lig2foam_vec_tmp_,
+        0.,
+        3 * lig2foam_->ncollected_
+    );
+    allocateArray
+    (
+        lig2foam_scl_tmp_,
+        0.,
+        lig2foam_->ncollected_
+    );
+}
+
+
+void twoWayOne2One::locateParticles()
+{
+    // get positions for locate
+    double** my_positions = static_cast<double**>(lmp->atom->x);
+    double*  my_flattened_positions = nullptr;
+    allocateArray(my_flattened_positions, 0., 3*lmp->atom->nlocal);
+    for (int atomi = 0; atomi < lmp->atom->nlocal; atomi++)
+    {
+        for (int coordi = 0; coordi < 3; coordi++)
+        {
+            my_flattened_positions[atomi*3+coordi] = my_positions[atomi][coordi];
+        }
+    }
+
+    double* collected_flattened_positions = nullptr;
+    allocateArray(collected_flattened_positions, 0., 3*lig2foam_->ncollected_);
+
+    lig2foam_->exchange(my_flattened_positions, collected_flattened_positions, 3);
+    destroy(my_flattened_positions);
+
+    double* my_prev_cell_ids = nullptr;
+    double* prev_cell_ids = nullptr;
+    if (cellIdComm_)
+    {
+        my_prev_cell_ids = my_prev_cell_ids_fix_->vector_atom;
+        allocateArray(prev_cell_ids, -1, lig2foam_->ncollected_);
+        lig2foam_->exchange(my_prev_cell_ids, prev_cell_ids);
+    }
+
+    if (lig2foam_mask_)
+    {
+        delete [] lig2foam_mask_;
+    }
+    lig2foam_mask_ = new bool[lig2foam_->ncollected_];
+
+    DynamicList<label> cellIds;
+    cellIds.setCapacity(lig2foam_->ncollected_);
+    label n_located(0);
+    label roundedCelli(-1);
+    const label nCells(particleCloud_.mesh().cells().size());
+    for (int atomi = 0; atomi < lig2foam_->ncollected_; atomi++)
+    {
+        const vector position = vector
+        (
+            collected_flattened_positions[3*atomi+0],
+            collected_flattened_positions[3*atomi+1],
+            collected_flattened_positions[3*atomi+2]
+        );
+        if (!thisFoamBox_.contains(position))
+        {
+            lig2foam_mask_[atomi] = false;
+            continue;
+        }
+        const label cellI = particleCloud_.locateM().findSingleCell
+        (
+            position,
+            cellIdComm_
+            ?   // don't know whether using round is efficient
+                (roundedCelli = round(prev_cell_ids[atomi])) < nCells
+                ?
+                roundedCelli
+                :
+                -1
+            :
+            -1
+        );
+
+        lig2foam_mask_[atomi] = false;
+        if (cellI >= 0) // in domain
+        {
+            lig2foam_mask_[atomi] = true;
+            n_located++;
+            cellIds.append(cellI);
+        }
+    }
+    if (cellIdComm_)
+    {
+        destroy(prev_cell_ids);
+    }
+
+    setNumberOfParticles(n_located);
+    particleCloud_.reAllocArrays();
+
+    reduce(n_located, sumOp<label>());
+    if (verbose_ || n_located != returnReduce(lmp->atom->nlocal, sumOp<label>()))
+    {
+        Warning << "Have located " << n_located
+                << " ouf of " << returnReduce(lmp->atom->nlocal, sumOp<label>())
+                << " particles in FOAM. "
+                << endl;
+    }
+
+    // copy positions/cellids/ids of located particles into arrays
+    allocateArray(lig2foam_ids_, 0, getNumberOfParticles());
+    int* collected_ids = nullptr;
+    allocateArray(collected_ids, 0, lig2foam_->ncollected_);
+    lig2foam_->exchange<int>(lmp->atom->tag, collected_ids);
+    extractCollected<int>(collected_ids, lig2foam_ids_);
+    destroy(collected_ids);
+
+    double* extracted_flattened_positions = new double[getNumberOfParticles()*3];
+    extractCollected<double>
+    (
+        collected_flattened_positions,
+        extracted_flattened_positions,
+        3
+    );
+    setPositions(getNumberOfParticles(), extracted_flattened_positions);
+    destroy(extracted_flattened_positions);
+    destroy(collected_flattened_positions);
+
+    setCellIDs(cellIds);
+}
+
+void twoWayOne2One::setupFoam2LigCommunication()
+{
+    int* src_procs = new int[thisFoamPartner_.size()];
+    for (int proci = 0; proci < thisFoamPartner_.size(); proci++)
+    {
+        src_procs[proci] = thisFoamPartner_[proci];
+    }
+
+    int* dst_procs = new int[thisLigPartner_.size()];
+    for (int proci = 0; proci < thisLigPartner_.size(); proci++)
+    {
+        dst_procs[proci] = thisLigPartner_[proci];
+    }
+
+    delete foam2lig_;
+    foam2lig_ = new One2One(comm_liggghts_);
+
+    foam2lig_->setup
+    (
+        thisFoamPartner_.size(),
+        src_procs,
+        thisLigPartner_.size(),
+        dst_procs,
+        getNumberOfParticles()
+    );
+    allocateArray
+    (
+        foam2lig_ids_,
+        0,
+        foam2lig_->ncollected_
+    );
+    foam2lig_->exchange<int>(lig2foam_ids_, foam2lig_ids_);
+
+    allocateArray
+    (
+        foam2lig_vec_tmp_,
+        0.,
+        3 * foam2lig_->ncollected_
+    );
+    allocateArray
+    (
+        foam2lig_scl_tmp_,
+        0.,
+        foam2lig_->ncollected_
+    );
+
+    if (cellIdComm_)
+    {
+        double** dbl_cell_ids = new double*[getNumberOfParticles()];
+        for (int atomi = 0; atomi < getNumberOfParticles(); atomi++)
+        {   // TEMPORARY: if this persists after 19.07.2018, call me.
+            dbl_cell_ids[atomi] = new double[1];
+            dbl_cell_ids[atomi][0] = particleCloud_.cellIDs()[atomi][0];
+        }
+        giveData("prev_cell_ids", "scalar-atom", dbl_cell_ids, "double");
+        delete [] dbl_cell_ids;
+    }
+}
+
+template <typename T>
+void twoWayOne2One::extractCollected(T**& src, T**& dst, int width) const
+{
+    int locali = 0;
+
+    for (int atomi = 0; atomi < lig2foam_->ncollected_; atomi++)
+    {
+        if (!lig2foam_mask_[atomi]) continue;
+
+        for (int coordi = 0; coordi < width; coordi++)
+        {
+            dst[locali][coordi] = src[atomi][coordi];
+        }
+        locali++;
+    }
+}
+
+template <typename T>
+void twoWayOne2One::extractCollected(T*& src, T*& dst, int width) const
+{
+    int locali = 0;
+
+    for (int atomi = 0; atomi < lig2foam_->ncollected_; atomi++)
+    {
+        if (!lig2foam_mask_[atomi]) continue;
+
+        for (int coordi = 0; coordi < width; coordi++)
+        {
+            dst[locali] = src[atomi*width+coordi];
+            locali++;
+        }
+    }
+}
+
+template <typename T>
+void twoWayOne2One::extractCollected(T*& src, T**& dst, int width) const
+{
+    int locali = 0;
+
+    for (int atomi = 0; atomi < lig2foam_->ncollected_; atomi++)
+    {
+        if (!lig2foam_mask_[atomi]) continue;
+
+        for (int coordi = 0; coordi < width; coordi++)
+        {
+            dst[locali][coordi] = src[atomi*width+coordi];
+        }
+        locali++;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,233 @@
+/*---------------------------------------------------------------------------*\
+    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
+    enhanced two way DEM-CFD coupling via MPI.
+
+    Compared to twoWayMPI, no Allreduces are used for communication.
+    Instead, a geometric map between FOAM and LIG domains is created and
+    subsequently used for communication.
+
+Class
+    twoWayOne2One
+
+SourceFiles
+    twoWayOne2One.C
+
+Contributing authors
+    Paul Kieckhefen (TUHH)  2018
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef twoWayOne2One_H
+#define twoWayOne2One_H
+
+#include "dataExchangeModel.H"
+#include "liggghtsCommandModel.H"
+#include "OFstream.H"
+#include <sys/stat.h>
+#include "pair.h"
+#include "force.h"
+#include "forceModel.H"
+#include "one2one.H"
+
+#include "meshSearch.H"
+
+//=================================//
+//LAMMPS/LIGGGHTS
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <mpi.h>
+#include <lammps.h>         // these are LAMMPS include files
+#include <input.h>
+#include <atom.h>
+#include <memory.h>
+#include <library.h>
+#include <library_cfd_coupling.h>
+#include <update.h>
+#include <comm.h>
+#include <fix.h>
+#include <fix_property_atom.h>
+//=================================//
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class noDrag Declaration
+\*---------------------------------------------------------------------------*/
+
+class twoWayOne2One
+:
+    public dataExchangeModel
+{
+private:
+
+  // private data
+    dictionary        propsDict_;
+
+    MPI_Comm  comm_liggghts_;
+
+
+    // LIG ranks from which to retrieve particle data
+    labelList thisLigPartner_;
+    labelList thisFoamPartner_;
+
+    One2One*  lig2foam_;
+    One2One*  foam2lig_;
+
+    bool*     lig2foam_mask_;
+
+    int*      lig2foam_ids_;
+    int*      foam2lig_ids_;
+
+    mutable double*   lig2foam_vec_tmp_;
+    mutable double*   lig2foam_scl_tmp_;
+
+    mutable double*   foam2lig_vec_tmp_;
+    mutable double*   foam2lig_scl_tmp_;
+
+    Switch    staticProcMap_;
+    Switch    cellIdComm_;
+    LAMMPS_NS::FixPropertyAtom* my_prev_cell_ids_fix_;
+
+    treeBoundBox       thisFoamBox_;
+
+    Switch            verbose_;
+
+
+  // private member functions
+
+    //- creates a geometric mapping between FOAM and LIG domains
+    void createProcMap();
+
+    //- create a One2One communicator which transfers from LIG to FOAM
+    void setupLig2FoamCommunication();
+
+    //- locates particles received from Lig
+    void locateParticles();
+
+    //- create a One2One communicator which transfers from FOAM to LIG
+    void setupFoam2LigCommunication();
+
+protected:
+    LAMMPS_NS::LAMMPS *lmp;
+
+public:
+
+    //- Runtime type information
+    TypeName("twoWayOne2One");
+
+
+    // Constructors
+
+        //- Construct from components
+        twoWayOne2One
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~twoWayOne2One();
+
+
+    // Member Functions
+        void getData
+        (
+            word name,
+            word type,
+            double ** const& field,
+            label step
+        ) const;
+
+        void getData
+        (
+            word name,
+            word type,
+            int ** const& field,
+            label step
+        ) const;
+
+        void giveData
+        (
+            word name,
+            word type,
+            double ** const& field,
+            const char* datatype
+        ) const;
+
+        //============
+        // double **
+        void allocateArray(double**&, double, int, int) const;
+        void allocateArray(double**&, double, int,const char* = "nparticles") const;
+        void inline destroy(double**,int=0) const;
+        //============
+        // int **
+        void allocateArray(int**&, int, int, int) const;
+        void allocateArray(int**&, int, int,const char* = "nparticles") const;
+        void inline destroy(int**,int=0) const;
+        //==============
+
+        //==============
+        // double *
+        void allocateArray(double*&, double, int) const;
+        void inline destroy(double*) const;
+        //==============
+        // int *
+        void allocateArray(int*&, int, int) const;
+        void inline destroy(int*) const;
+        //==============
+
+        bool couple(int);
+
+        //- extractCollected takes the collected data from Lig
+        //  present in this Foam domain and applies the mask.
+        //  the width parameter can be used for reshaping.
+        template <typename T>
+        void extractCollected(T**&, T**&, int width=1) const;
+        template <typename T>
+        void extractCollected(T*&, T*&, int width=1) const;
+        template <typename T>
+        void extractCollected(T*&, T**&, int width=1) const;
+
+        scalar getCG() const { return lmp->force->cg(); }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -74,7 +74,7 @@ energyModel::~energyModel()
 
 scalar energyModel::Cp() const
 {
-    return Cp_; 
+    return Cp_;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/energyModel/energyModel/energyModel.H

@@ -43,13 +43,13 @@ protected:
         const dictionary& dict_;
 
         cfdemCloudEnergy& particleCloud_;
-	
-	IOdictionary transportProperties_;
-	
-	scalar kf0_;                              // fluid thermal conductivity [W/(m*K)]
 
-        scalar Cp_;                                  // specific heat capacity [W*s/(kg*K)]
-       
+        IOdictionary transportProperties_;
+
+        scalar kf0_;                                // fluid thermal conductivity [W/(m*K)]
+
+        scalar Cp_;                                 // specific heat capacity [W*s/(kg*K)]
+
 public:
 
     //- Runtime type information
@@ -76,7 +76,7 @@ public:
         energyModel
         (
             const dictionary& dict,
-            cfdemCloudEnergy& sm            
+            cfdemCloudEnergy& sm
         );
 
 
@@ -98,16 +98,16 @@ public:
     // Member Functions
 
         virtual void addEnergyContribution(volScalarField&) const = 0;
-	
-	virtual void addEnergyCoefficient(volScalarField&) const = 0;
-	
-	virtual void calcEnergyContribution() = 0;
-	
-	virtual void postFlow() {}
 
-	virtual void solve() {}
-	
-	scalar Cp() const;
+        virtual void addEnergyCoefficient(volScalarField&) const = 0;
+
+        virtual void calcEnergyContribution() = 0;
+
+        virtual void postFlow() {}
+
+        virtual void solve() {}
+
+        scalar Cp() const;
 
 
 };

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

@@ -21,6 +21,7 @@ License
 #include "error.H"
 #include "heatTransferGunn.H"
 #include "addToRunTimeSelectionTable.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -43,6 +44,7 @@ heatTransferGunn::heatTransferGunn
 :
     energyModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
     expNusselt_(propsDict_.lookupOrDefault<bool>("expNusselt",false)),
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
@@ -141,7 +143,9 @@ heatTransferGunn::heatTransferGunn
     partHeatFlux_(NULL),
     partHeatFluxCoeff_(NULL),
     partRe_(NULL),
-    partNu_(NULL)
+    partNu_(NULL),
+    scaleDia_(1.),
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0)))
 {
     allocateMyArrays();
 
@@ -181,6 +185,14 @@ heatTransferGunn::heatTransferGunn
           FatalError <<"Cannot read and create NuField at the same time!\n" << abort(FatalError);
         }
     }
+
+    if (propsDict_.found("scale") && typeCG_.size()==1)
+    {
+        // if "scale" is specified and there's only one single type, use "scale"
+        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        typeCG_[0] = scaleDia_;
+    }
+    else if (typeCG_.size()>1) multiTypes_ = true;
 }
 
 
@@ -255,6 +267,15 @@ void heatTransferGunn::calcEnergyContribution()
        const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
     #endif
 
+    if (typeCG_.size()>1 || typeCG_[0] > 1)
+    {
+        Info << "heatTransferGunn using scale = " << typeCG_ << endl;
+    }
+    else if (particleCloud_.cg() > 1)
+    {
+        scaleDia_=particleCloud_.cg();
+        Info << "heatTransferGunn using scale from liggghts cg = " << scaleDia_ << endl;
+    }
 
     // calc La based heat flux
     scalar voidfraction(1);
@@ -263,6 +284,8 @@ void heatTransferGunn::calcEnergyContribution()
     label cellI=0;
     vector Us(0,0,0);
     scalar ds(0);
+    scalar ds_scaled(0);
+    scalar scaleDia3 = typeCG_[0]*typeCG_[0]*typeCG_[0];
     scalar muf(0);
     scalar magUr(0);
     scalar Rep(0);
@@ -270,6 +293,8 @@ void heatTransferGunn::calcEnergyContribution()
     scalar Nup(0);
     scalar Tsum(0.0);
 
+    scalar cg = typeCG_[0];
+    label partType = 1;
 
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
@@ -297,22 +322,30 @@ void heatTransferGunn::calcEnergyContribution()
                 if (voidfraction < 0.01)
                     voidfraction = 0.01;
 
+                if (multiTypes_)
+                {
+                    partType = particleCloud_.particleType(index);
+                    cg = typeCG_[partType - 1];
+                    scaleDia3 = cg*cg*cg;
+                }
+
                 // calc relative velocity
                 Us = particleCloud_.velocity(index);
                 magUr = mag(Ufluid - Us);
                 ds = 2.*particleCloud_.radius(index);
+                ds_scaled = ds/cg;
                 muf = mufField[cellI];
-                Rep = ds * magUr * voidfraction * rho_[cellI]/ muf;
+                Rep = ds_scaled * magUr * voidfraction * rho_[cellI]/ muf;
                 Pr = max(SMALL, Cp_ * muf / kf0_);
 
                 Nup = Nusselt(voidfraction, Rep, Pr);
 
                 Tsum += partTemp_[index][0];
-                scalar h = kf0_ * Nup / ds;
-                scalar As = ds * ds * M_PI; // surface area of sphere
+                scalar h = kf0_ * Nup / ds_scaled;
+                scalar As = ds_scaled * ds_scaled * M_PI; // surface area of sphere
 
                 // calc convective heat flux [W]
-                heatFlux(index, h, As, Tfluid);
+                heatFlux(index, h, As, Tfluid, scaleDia3);
 
                 if(verbose_)
                 {
@@ -434,9 +467,16 @@ scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) con
                         Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33);
 }
 
-void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid)
+void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid, scalar cg3)
 {
-    scalar hAs = h * As;
+    scalar hAs = h * As * cg3;
+
+    if (particleCloud_.getParticleEffVolFactors())
+    {
+        scalar effVolFac = particleCloud_.particleEffVolFactor(index);
+        hAs *= effVolFac;
+    }
+
     partHeatFlux_[index][0] = - hAs * partTemp_[index][0];
     if(!implicit_)
     {

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

@@ -28,6 +28,7 @@ License
 #include "fvCFD.H"
 #include "cfdemCloudEnergy.H"
 #include "energyModel.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -45,6 +46,8 @@ protected:
 
     dictionary propsDict_;
 
+    bool multiTypes_;
+
     bool expNusselt_;
 
     bool interpolation_;
@@ -79,13 +82,13 @@ protected:
 
     word tempFieldName_;
 
-    const volScalarField& tempField_;            // ref to temperature field
+    const volScalarField& tempField_;           // ref to temperature field
 
     word voidfractionFieldName_;
 
-    const volScalarField& voidfraction_;    // ref to voidfraction field
+    const volScalarField& voidfraction_;        // ref to voidfraction field
 
-    scalar maxSource_;                           // max (limited) value of src field
+    scalar maxSource_;                          // max (limited) value of src field
 
     word velFieldName_;
 
@@ -109,6 +112,10 @@ protected:
 
     mutable double **partNu_;
 
+    mutable scalar scaleDia_;
+
+    scalarList typeCG_;
+
     void allocateMyArrays() const;
 
     void partTempField();
@@ -117,7 +124,7 @@ protected:
 
     virtual void giveData();
 
-    virtual void heatFlux(label, scalar, scalar, scalar);
+    virtual void heatFlux(label, scalar, scalar, scalar, scalar cg3 = 1.0);
 
 public:
 
@@ -143,11 +150,11 @@ public:
 
         void addEnergyContribution(volScalarField&) const;
 
-	void addEnergyCoefficient(volScalarField&) const;
+        void addEnergyCoefficient(volScalarField&) const;
 
-	void calcEnergyContribution();
+        void calcEnergyContribution();
 
-	void postFlow();
+        void postFlow();
 
         scalar aveTpart() const;
 };

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

@@ -45,7 +45,7 @@ heatTransferGunnPartField::heatTransferGunnPartField
 :
     heatTransferGunn(dict,sm),
     partCpField_
-    (   
+    (
         IOobject
         (
             "partCp",
@@ -82,15 +82,15 @@ heatTransferGunnPartField::heatTransferGunnPartField
     {
         FatalError << "heatTransferGunnPartField: provide list of specific heat capacities." << abort(FatalError);
     }
-    
+
     if (propsDict_.found("pTMax"))
     {
-        pTMax_.value()=scalar(readScalar(propsDict_.lookup("pTMax"))); 
+        pTMax_.value()=scalar(readScalar(propsDict_.lookup("pTMax")));
     }
-    
+
     if (propsDict_.found("pTMin"))
     {
-        pTMin_.value()=scalar(readScalar(propsDict_.lookup("pTMin"))); 
+        pTMin_.value()=scalar(readScalar(propsDict_.lookup("pTMin")));
     }
 
     partTempField_.writeOpt() = IOobject::AUTO_WRITE;

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H

@@ -15,7 +15,7 @@ License
     along with this code.  If not, see <http://www.gnu.org/licenses/>.
 
     Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
-    
+
     Description
     Correlation for Nusselt number according to
     Gunn, D. J. International Journal of Heat and Mass Transfer 21.4 (1978)
@@ -30,6 +30,7 @@ License
 #include "heatTransferGunn.H"
 #include "fvOptions.H"
 #include "scalarList.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -88,9 +89,9 @@ public:
     // Member Functions
         void addEnergyContribution(volScalarField&) const;
 
-	void calcEnergyContribution();
+        void calcEnergyContribution();
 
-	void postFlow();
+        void postFlow();
 
         void solve();
 

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

@@ -22,6 +22,7 @@ License
 #include "reactionHeat.H"
 #include "addToRunTimeSelectionTable.H"
 #include "dataExchangeModel.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H

@@ -24,6 +24,7 @@ License
 #include "fvCFD.H"
 #include "cfdemCloudEnergy.H"
 #include "energyModel.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -66,24 +67,24 @@ public:
 
     // Constructors
 
-    //- Construct from components
-    reactionHeat
-    (
-        const dictionary& dict,
-        cfdemCloudEnergy& sm
-    );
+        //- Construct from components
+        reactionHeat
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
 
     // Destructor
 
-    virtual ~reactionHeat();
+        virtual ~reactionHeat();
 
     // Member Functions
 
-    void addEnergyContribution(volScalarField&) const;
+        void addEnergyContribution(volScalarField&) const;
 
-    void addEnergyCoefficient(volScalarField&) const {}
+        void addEnergyCoefficient(volScalarField&) const {}
 
-    void calcEnergyContribution();
+        void calcEnergyContribution();
 };
 
 

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

@@ -0,0 +1,247 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "gradPForceSmooth.H"
+#include "addToRunTimeSelectionTable.H"
+#include "smoothingModel.H"
+#include "fvCFD.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(gradPForceSmooth, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    gradPForceSmooth,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+gradPForceSmooth::gradPForceSmooth
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    pFieldName_(propsDict_.lookup("pFieldName")),
+    p_(sm.mesh().lookupObject<volScalarField> (pFieldName_)),
+    velocityFieldName_(propsDict_.lookup("velocityFieldName")),
+    U_(sm.mesh().lookupObject<volVectorField> (velocityFieldName_)),
+    useRho_(false),
+    useU_(false),
+    addedMassCoeff_(0.0),
+    smoothingModel_
+    (
+        smoothingModel::New
+        (
+            propsDict_,
+            sm
+        )
+    ),
+    pSmooth_
+    (
+        IOobject
+        (
+            "pSmooth",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh()
+    )
+{
+    // init force sub model
+    setForceSubModels(propsDict_);
+
+    // define switches which can be read from dict
+    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch
+    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+
+    // read those switches defined above, if provided in dict
+    forceSubM(0).readSwitches();
+
+    if (modelType_ == "B")
+    {
+        FatalError <<"using  model gradPForceSmooth with model type B is not valid\n" << abort(FatalError);
+    }else if (modelType_ == "Bfull")
+    {
+        if(forceSubM(0).switches()[1])
+        {
+            Info << "Using treatForceDEM false!" << endl;
+            forceSubM(0).setSwitches(1,false); // treatForceDEM = false
+        }
+    }else // modelType_=="A"
+    {
+        if(!forceSubM(0).switches()[1])
+        {
+            Info << "Using treatForceDEM true!" << endl;
+            forceSubM(0).setSwitches(1,true); // treatForceDEM = true
+        }
+    }
+
+    if (propsDict_.found("useU")) useU_=true;
+    if (propsDict_.found("useAddedMass")) 
+    {
+        addedMassCoeff_ =  readScalar(propsDict_.lookup("useAddedMass"));
+        Info << "gradP will also include added mass with coefficient: " << addedMassCoeff_ << endl;
+        Info << "WARNING: use fix nve/sphere/addedMass in LIGGGHTS input script to correctly account for added mass effects!" << endl;
+    }
+
+    if(p_.dimensions()==dimensionSet(0,2,-2,0,0))
+        useRho_ = true;
+
+    particleCloud_.checkCG(true);
+
+    particleCloud_.probeM().initialize(typeName, "gradP.logDat");
+    particleCloud_.probeM().vectorFields_.append("gradPForceSmooth"); //first entry must the be the force
+    particleCloud_.probeM().scalarFields_.append("Vs");
+    particleCloud_.probeM().scalarFields_.append("rho");
+    particleCloud_.probeM().writeHeader();
+    
+    pSmooth_ = p_;
+   }
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+gradPForceSmooth::~gradPForceSmooth()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void gradPForceSmooth::setForce() const
+{
+    volVectorField gradPField = fvc::grad(p_);
+    if(pFieldName_ == "p_rgh")
+    {
+        const volScalarField& rho_ = particleCloud_.mesh().lookupObject<volScalarField>("rho");
+        const volScalarField& gh_ = particleCloud_.mesh().lookupObject<volScalarField>("gh");
+        
+        //Smooth p_rgh, easier to handle boundaries
+        smoothingM().smoothen(pSmooth_);
+        
+        //Superpose hydrostatic pressure
+        volScalarField pFull = pSmooth_ + rho_*gh_;
+        
+        gradPField = fvc::grad(pFull);
+        
+    }else{
+        
+        smoothingM().smoothen(pSmooth_);
+        gradPField = fvc::grad(pSmooth_);
+        
+    }
+    
+    /*if (useU_)
+    {
+        // const volScalarField& voidfraction_ = particleCloud_.mesh().lookupObject<volScalarField> ("voidfraction");
+        volScalarField U2 = U_&U_;// *voidfraction_*voidfraction_;
+        if (useRho_)
+            gradPField = fvc::grad(0.5*U2);
+        else
+            gradPField = fvc::grad(0.5*forceSubM(0).rhoField()*U2);
+    }*/
+    vector gradP;
+    scalar Vs;
+    scalar rho;
+    vector position;
+    vector force;
+    label cellI;
+
+    interpolationCellPoint<vector> gradPInterpolator_(gradPField);
+    #include "setupProbeModel.H"
+
+    for(int index = 0;index <  particleCloud_.numberOfParticles(); index++)
+    {
+        //if(mask[index][0])
+        //{
+            force=vector(0,0,0);
+            cellI = particleCloud_.cellIDs()[index][0];
+
+            if (cellI > -1) // particle Found
+            {
+                position = particleCloud_.position(index);
+                
+                if(forceSubM(0).interpolation()) // use intepolated values for alpha (normally off!!!)
+                {
+                    gradP = gradPInterpolator_.interpolate(position,cellI);
+                }else
+                {
+                    gradP = gradPField[cellI];
+                }
+
+                Vs = particleCloud_.particleVolume(index);
+                rho = forceSubM(0).rhoField()[cellI];
+
+                // calc particle's pressure gradient force
+                if (useRho_)
+                    force = -Vs*gradP*rho*(1.0+addedMassCoeff_);
+                else
+                    force = -Vs*gradP*(1.0+addedMassCoeff_);
+
+                if(forceSubM(0).verbose() && index >=0 && index <2)
+                {
+                    Info << "index = " << index << endl;
+                    Info << "gradP = " << gradP << endl;
+                    Info << "force = " << force << endl;
+                }
+
+                //Set value fields and write the probe
+                if(probeIt_)
+                {
+                    #include "setupProbeModelfields.H"
+                    vValues.append(force);           //first entry must the be the force
+                    sValues.append(Vs);
+                    sValues.append(rho);
+                    particleCloud_.probeM().writeProbe(index, sValues, vValues);
+                }
+            }
+
+            // write particle based data to global array
+            forceSubM(0).partToArray(index,force,vector::zero);
+
+        //}
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Class
+    gradPForceSmooth
+
+Description
+    This code is an extention of the gradPForce class by allowing the pressure-
+    field to be smoothened prior to the force calculation (without altering the
+    original field).
+
+SourceFiles
+    gradPForceSmooth.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef gradPForceSmooth_H
+#define gradPForceSmooth_H
+
+#include "forceModel.H"
+#include "interpolationCellPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                           Class gradPForceSmooth Declaration
+\*---------------------------------------------------------------------------*/
+
+class gradPForceSmooth
+:
+    public forceModel
+{
+private:
+    dictionary propsDict_;
+
+    word pFieldName_;
+
+    const volScalarField& p_;
+
+    word velocityFieldName_;
+
+    const volVectorField& U_;
+
+    bool useRho_;
+
+    bool useU_;          // if false: substitution p=0.5*rho*U^2
+
+    mutable double addedMassCoeff_; //added mass coefficient
+    
+    autoPtr<smoothingModel> smoothingModel_;
+    
+    mutable volScalarField pSmooth_;    
+
+public:
+
+    //- Runtime type information
+    TypeName("gradPForceSmooth");
+
+
+    // Constructors
+
+        //- Construct from components
+        gradPForceSmooth
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~gradPForceSmooth();
+
+
+    // Member Functions
+        void setForce() const;
+        inline const smoothingModel& smoothingM() const
+        {
+            return smoothingModel_;
+        }        
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code. If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "surfaceTensionForce.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(surfaceTensionForce, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    surfaceTensionForce,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+surfaceTensionForce::surfaceTensionForce
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    stfFieldName_(propsDict_.lookupOrDefault<word>("stfFieldName", "surfaceTensionForce")),
+    stf_(sm.mesh().lookupObject<surfaceScalarField> (stfFieldName_))
+{
+
+
+    // init force sub model
+    setForceSubModels(propsDict_);
+
+    // define switches which can be read from dict
+    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
+
+    // read those switches defined above, if provided in dict
+    forceSubM(0).readSwitches();
+
+    Info << "check if interpolation really works - use directly interpolationCellPoint<vector> ???" << endl;
+    particleCloud_.checkCG(false);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+surfaceTensionForce::~surfaceTensionForce()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void surfaceTensionForce::setForce() const
+{
+    volVectorField reconstructedStf = fvc::reconstruct(stf_*particleCloud_.mesh().magSf());
+    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    {
+        //if(mask[index][0])
+        //{
+            // definition of spherical particle
+            label cellI = particleCloud_.cellIDs()[index][0];
+            scalar Vp = particleCloud_.particleVolume(index);
+            if(cellI >-1.0) // particle found on proc domain
+            {
+                vector surfaceTensionForcep = Foam::vector(0,0,0);
+                
+                surfaceTensionForcep = Vp * reconstructedStf[cellI];
+
+               // write particle based data to global array
+               forceSubM(0).partToArray(index,surfaceTensionForcep,vector::zero);
+
+            } // end if particle found on proc domain
+        //}// end if in mask
+    }// end loop particles
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code. If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Description
+    Applies a pre-calculated surface tension force to the particles.
+
+Class
+    surfaceTensionForce
+
+SourceFiles
+    surfaceTensionForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceTensionForce_H
+#define surfaceTensionForce_H
+
+#include "forceModel.H"
+#include "interpolationCellPoint.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class surfaceTensionForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfaceTensionForce
+:
+    public forceModel
+{
+private:
+    dictionary propsDict_;
+
+    word stfFieldName_;
+    
+    const surfaceScalarField& stf_;
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaceTensionForce");
+
+
+    // Constructors
+
+        //- Construct from components
+        surfaceTensionForce
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~surfaceTensionForce();
+
+
+    // Member Functions
+        void setForce() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -95,7 +95,16 @@ label engineSearchMany2Many::intersection
     return face;
 }
 
-
+label engineSearchMany2Many::findCell
+(
+    double** const& mask,
+    double**& positions,
+    int**& cellIDs,
+    int size
+) const
+{
+    return 1; // locate is provided by many2may / one2one
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -87,6 +87,13 @@ public:
                 const point& pEnd
         ) const;
 
+        label findCell
+        (
+            double** const& mask,
+            double**& positions,
+            int**& cellIDs,
+            int size
+        ) const;
 };
 
 

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

@@ -19,6 +19,7 @@ License
 #include "expParticleForces.H"
 #include "mathExtra.H"
 #include "addToRunTimeSelectionTable.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

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

@@ -26,6 +26,7 @@ SourceFiles
 #define expParticleForces_H
 
 #include "otherForceModel.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

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

@@ -19,6 +19,7 @@ License
 #include "gravity.H"
 #include "mathExtra.H"
 #include "addToRunTimeSelectionTable.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

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

@@ -26,6 +26,7 @@ SourceFiles
 #define gravity_H
 
 #include "otherForceModel.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

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

@@ -18,6 +18,7 @@ License
 #include "error.H"
 #include "otherForceModel.H"
 #include "mathExtra.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

src/lagrangian/cfdemParticle/subModels/otherForceModel/otherForceModel/otherForceModel.H

@@ -28,6 +28,7 @@ SourceFiles
 #include "fvCFD.H"
 #include "cfdemCloud.H"
 #include "probeModel.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

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

@@ -19,6 +19,7 @@ License
 #include "weightSecondaryPhase.H"
 #include "mathExtra.H"
 #include "addToRunTimeSelectionTable.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -82,11 +83,11 @@ tmp<volVectorField> weightSecondaryPhase::exportForceField()
             )
         )
     );
-    
+
     volVectorField& source = tsource.ref();
-    
+
     source = rho_ * alpha_ * g_;
-    
+
     return tsource;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -26,6 +26,7 @@ SourceFiles
 #define weightSecondaryPhase_H
 
 #include "otherForceModel.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -48,12 +49,12 @@ protected:
         word volfracFieldName_;
 
         const volScalarField& alpha_;
-    
+
         dimensionedScalar rho_;
 
         const dimensionedVector g_;
 
-    
+
 public:
 
     //- Runtime type information
@@ -74,7 +75,6 @@ public:
         virtual ~weightSecondaryPhase();
 
 
-
     // Member Functions
 
 

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

@@ -6,7 +6,7 @@
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
                                 Copyright 2012-     DCS Computing GmbH, Linz
-                                Copyright (C) 2013-     Graz University of  
+                                Copyright (C) 2013-     Graz University of
                                                         Technology, IPPT
 -------------------------------------------------------------------------------
 License
@@ -91,7 +91,7 @@ public:
         bool doSmoothing() const;
 
         //void dSmoothing(volScalarField&) const;
-        
+
         void smoothen(volScalarField&) const;
 
         void smoothen(volVectorField&) const;

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

@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "temporalSmoothing.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(temporalSmoothing, 0);
+
+addToRunTimeSelectionTable
+(
+    smoothingModel,
+    temporalSmoothing,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+temporalSmoothing::temporalSmoothing
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    smoothingModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    lowerLimit_(readScalar(propsDict_.lookup("lowerLimit"))),
+    upperLimit_(readScalar(propsDict_.lookup("upperLimit"))),
+    verbose_(false),
+    refFieldName_(propsDict_.lookup("refField")),
+    gamma_(readScalar(propsDict_.lookup("gamma")))
+{
+
+    if(propsDict_.found("verbose"))
+        verbose_ = true;
+
+    checkFields(sSmoothField_);
+    checkFields(vSmoothField_);
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+temporalSmoothing::~temporalSmoothing()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+bool temporalSmoothing::doSmoothing() const
+{
+    return true;
+}
+
+
+void Foam::temporalSmoothing::smoothen(volScalarField& fieldSrc) const
+{
+    // Create scalar smooth field from virgin scalar smooth field template
+    volScalarField sSmoothField = sSmoothField_;
+
+    sSmoothField.dimensions().reset(fieldSrc.dimensions());
+    sSmoothField.ref()=fieldSrc.internalField();
+    sSmoothField.correctBoundaryConditions();
+    sSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
+    sSmoothField.oldTime()=fieldSrc;
+    sSmoothField.oldTime().correctBoundaryConditions();
+
+    volScalarField refField = particleCloud_.mesh().lookupObject<volScalarField>(refFieldName_);
+
+    // do smoothing
+    dimensionedScalar deltaT = sSmoothField.mesh().time().deltaT();
+    solve
+    (
+        fvm::ddt(sSmoothField)
+        -
+        gamma_/deltaT * (refField - fvm::Sp(1.0,sSmoothField))
+    );
+
+    // bound sSmoothField_
+    forAll(sSmoothField,cellI)
+    {
+        sSmoothField[cellI]=max(lowerLimit_,min(upperLimit_,sSmoothField[cellI]));
+    }
+
+    // get data from working sSmoothField - will copy only values at new time
+    fieldSrc=sSmoothField;
+    fieldSrc.correctBoundaryConditions();
+
+    if(verbose_)
+    {
+        Info << "min/max(fieldoldTime) (unsmoothed): " << min(sSmoothField.oldTime()) << tab << max(sSmoothField.oldTime()) << endl;
+        Info << "min/max(fieldSrc): " << min(fieldSrc) << tab << max(fieldSrc) << endl;
+        Info << "min/max(fieldSrc.oldTime): " << min(fieldSrc.oldTime()) << tab << max(fieldSrc.oldTime()) << endl;
+    }
+
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::temporalSmoothing::smoothen(volVectorField& fieldSrc) const
+{
+    // Create scalar smooth field from virgin scalar smooth field template
+    volVectorField vSmoothField = vSmoothField_;
+
+    vSmoothField.dimensions().reset(fieldSrc.dimensions());
+    vSmoothField.ref()=fieldSrc.internalField();
+    vSmoothField.correctBoundaryConditions();
+    vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
+    vSmoothField.oldTime()=fieldSrc;
+    vSmoothField.oldTime().correctBoundaryConditions();
+
+    volVectorField refField = particleCloud_.mesh().lookupObject<volVectorField>(refFieldName_);
+
+    dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
+    solve
+    (
+        fvm::ddt(vSmoothField)
+        -
+        gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
+    );
+
+    // get data from working vSmoothField
+    fieldSrc=vSmoothField;
+    fieldSrc.correctBoundaryConditions();
+
+    if(verbose_)
+    {
+        Info << "min/max(fieldoldTime) (unsmoothed): " << min(vSmoothField.oldTime()) << tab << max(vSmoothField.oldTime()) << endl;
+        Info << "min/max(fieldSrc): " << min(fieldSrc) << tab << max(fieldSrc) << endl;
+        Info << "min/max(fieldSrc.oldTime): " << min(fieldSrc.oldTime()) << tab << max(fieldSrc.oldTime()) << endl;
+    }
+
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::temporalSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const
+{
+    FatalError << "Smoothen reference field is not implemented for this smoothing model!" << abort(FatalError);
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/smoothingModel/temporalSmoothing/temporalSmoothing.H

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Class
+    temporalSmoothing
+
+Description
+    Smoothens a field using a temporal relaxation approach with a specified
+    reference field.
+
+SourceFiles
+    temporalSmoothing.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef temporalSmoothing_H
+#define temporalSmoothing_H
+
+#include "smoothingModel.H"
+//#include "wallFvPatch.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class temporalSmoothing Declaration
+\*---------------------------------------------------------------------------*/
+
+class temporalSmoothing
+:
+    public smoothingModel
+{
+
+private:
+
+    dictionary propsDict_;
+    scalar lowerLimit_;
+    scalar upperLimit_;
+    bool verbose_;
+    word refFieldName_;
+    scalar gamma_;
+
+public:
+
+    //- Runtime type information
+    TypeName("temporalSmoothing");
+
+
+    // Constructors
+
+        //- Construct from components
+        temporalSmoothing
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~temporalSmoothing();
+
+
+    // Member Functions
+        bool doSmoothing() const;
+
+        //void dSmoothing(volScalarField&) const;
+
+        void smoothen(volScalarField&) const;
+
+        void smoothen(volVectorField&) const;
+
+        void smoothenReferenceField(volVectorField&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -48,24 +48,24 @@ class SyamlalThermCond
 :
     public thermCondModel
 {
-  
+
 private:
-  
+
     dictionary propsDict_;
-  
+
     word voidfractionFieldName_;
-  
+
     const volScalarField& voidfraction_;
-    
+
     word rhoFieldName_;
-    
+
     const volScalarField& rho_;
-    
+
     word wallQFactorName_;
-    
+
     // ratio of half-cell-size and near-wall film
     mutable volScalarField wallQFactor_;
-    
+
     bool hasWallQFactor_;
 
 public:
@@ -89,12 +89,11 @@ public:
 
 
     // Member Functions
-	
-	tmp<volScalarField> thermCond() const;
-	
-	tmp<volScalarField> thermDiff() const;
-       
-       
+
+        tmp<volScalarField> thermCond() const;
+
+        tmp<volScalarField> thermDiff() const;
+
 };
 
 

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

@@ -62,7 +62,7 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
         ),
         sm.mesh(),
         dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
-	"zeroGradient"
+        "zeroGradient"
     ),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),

src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H

@@ -49,26 +49,26 @@ class ZehnerSchluenderThermCond
 :
     public thermCondModel
 {
-  
+
 private:
-  
+
     dictionary propsDict_;
 
     mutable volScalarField partKsField_;
-  
+
     word voidfractionFieldName_;
-  
+
     const volScalarField& voidfraction_;
-    
+
     scalarList typeKs_;
 
     mutable double **partKs_;
-    
+
     word wallQFactorName_;
-    
+
     // ratio of half-cell-size and near-wall film
     mutable volScalarField wallQFactor_;
-    
+
     bool hasWallQFactor_;
 
     void allocateMyArrays() const;
@@ -96,12 +96,11 @@ public:
 
 
     // Member Functions
-	
-	tmp<volScalarField> thermCond() const;
-	
-	tmp<volScalarField> thermDiff() const;
-       
-       
+
+        tmp<volScalarField> thermCond() const;
+
+        tmp<volScalarField> thermDiff() const;
+
 };
 
 

src/lagrangian/cfdemParticle/subModels/thermCondModel/noTherm/noThermCond.C

@@ -22,6 +22,7 @@ License
 
 #include "noThermCond.H"
 #include "addToRunTimeSelectionTable.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

src/lagrangian/cfdemParticle/subModels/thermCondModel/thermCondModel/thermCondModel.H

@@ -33,6 +33,7 @@ SourceFiles
 
 #include "fvCFD.H"
 #include "cfdemCloud.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -51,13 +52,13 @@ protected:
         const dictionary& dict_;
 
         cfdemCloud& particleCloud_;
-	
-	IOdictionary transportProperties_;
-	
-	dimensionedScalar kf0_;
-	
-	dimensionedScalar Cp_;
-	
+
+        IOdictionary transportProperties_;
+
+        dimensionedScalar kf0_;
+
+        dimensionedScalar Cp_;
+
 
 public:
 
@@ -104,11 +105,11 @@ public:
 
 
     // Member Functions
-	
-	virtual tmp<volScalarField> thermCond() const = 0;
-	
-	virtual tmp<volScalarField> thermDiff() const = 0;
-       
+
+        virtual tmp<volScalarField> thermCond() const = 0;
+
+        virtual tmp<volScalarField> thermDiff() const = 0;
+
 };
 
 

tutorials/cfdemSolverMultiphase/damBreak/Allclean.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+
+#- clean up case
+echo "deleting data at: $casePath :\n"
+source $WM_PROJECT_DIR/bin/tools/CleanFunctions
+cd $casePath/CFD
+cleanCase
+rm -r $casePath/CFD/clockData
+rm $casePath/DEM/post/*.*
+touch $casePath/DEM/post/.gitignore
+rm -r $casePath/CFD/0
+rm $casePath/log*
+
+echo "Remove restart file?"
+echo "Enter: yes, Ctrl + C: no"
+read
+
+rm $casePath/DEM/post/restart/*.*
+rm $casePath/DEM/post/restart/liggghts.restartCFDEM*

tutorials/cfdemSolverMultiphase/damBreak/Allrun.sh

@@ -0,0 +1,35 @@
+#!/bin/bash
+
+#===================================================================#
+# Allrun script for cfdemSolverMultiphase
+#===================================================================#
+
+#- define variables
+postProcessing=false #true
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+
+# check if mesh was built
+if [ -f "$casePath/CFD/constant/polyMesh/points" ]; then
+    echo "mesh was built before - using old mesh"
+else
+    echo "mesh needs to be built"
+    cd $casePath/CFD
+    blockMesh
+fi
+
+cd $casePath/CFD
+cp -r 0.org 0
+setFields
+
+if [ -f "$casePath/DEM/post/restart/liggghts.restart" ];  then
+    echo "LIGGGHTS init was run before - using existing restart file"
+else
+    #- run DEM in new terminal
+    $casePath/parDEMrun.sh
+fi
+
+bash $casePath/parCFDDEMrun.sh
+
+if [ "$postProcessing" = true ]; then
+    bash $casePath/postrun.sh
+fi