Merge pull request #146 from ParticulateFlow/release

Release 23.02

.gitignore

@@ -8,7 +8,7 @@ log.*
 *.swp
 *.swo
 
-**/linux*Gcc*/
+**/linux*cc*/
 **/.vscode
 
 lnInclude

applications/solvers/cfdemSolverIB/Make/options

@@ -9,12 +9,12 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -Wno-deprecated-copy
 
 EXE_LIBS = \
@@ -27,6 +27,7 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lfvOptions \
+    -lsampling \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -6,7 +6,8 @@
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright (C) 1991-2009 OpenCFD Ltd.
                                 Copyright (C) 2009-2012 JKU, Linz
-                                Copyright (C) 2012-     DCS Computing GmbH,Linz
+                                Copyright (C) 2012-2015 DCS Computing GmbH,Linz
+                                Copyright (C) 2015-     JKU, Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -29,11 +30,14 @@ Application
 
 Description
     Transient solver for incompressible flow.
-    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, 
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
     where additional functionality for CFD-DEM coupling using immersed body
     (fictitious domain) method is added.
 Contributions
     Alice Hager
+    Daniel Queteschiner
+    Thomas Lichtenegger
+    Achuth N. Balachandran Nair
 \*---------------------------------------------------------------------------*/
 
 
@@ -53,23 +57,21 @@ Contributions
 
 #include "cellSet.H"
 
+#include "fvOptions.H"  // added the fvOptions library
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+
     #include "setRootCase.H"
-
     #include "createTime.H"
-
     #include "createDynamicFvMesh.H"
-
     #include "createControl.H"
-
     #include "createTimeControls.H"
-
     #include "createFields.H"
-
     #include "initContinuityErrs.H"
+    #include "createFvOptions.H"
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
@@ -93,24 +95,31 @@ int main(int argc, char *argv[])
 
         // do particle stuff
         Info << "- evolve()" << endl;
-        particleCloud.evolve();
+        particleCloud.evolve(Us);
 
         // Pressure-velocity PISO corrector
         {
+            MRF.correctBoundaryVelocity(U);
+
             // Momentum predictor
 
             fvVectorMatrix UEqn
             (
-                fvm::ddt(voidfraction,U)
+                fvm::ddt(voidfraction,U) + MRF.DDt(U)
               + fvm::div(phi, U)
               + turbulence->divDevReff(U)
+             ==
+                fvOptions(U)
             );
 
             UEqn.relax();
 
+            fvOptions.constrain(UEqn);
+
             if (piso.momentumPredictor())
             {
                 solve(UEqn == -fvc::grad(p));
+                fvOptions.correct(U);
             }
 
             // --- PISO loop
@@ -126,6 +135,7 @@ int main(int argc, char *argv[])
 
                 adjustPhi(phi, U, p);
 
+
                 while (piso.correctNonOrthogonal())
                 {
                     // Pressure corrector
@@ -152,12 +162,15 @@ int main(int argc, char *argv[])
             }
         }
 
+        laminarTransport.correct();
         turbulence->correct();
 
         Info << "particleCloud.calcVelocityCorrection() " << endl;
         volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
         particleCloud.calcVelocityCorrection(p,U,phiIB,voidfractionNext);
 
+        fvOptions.correct(U);
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/cfdemSolverIB/createFields.H

@@ -26,21 +26,6 @@
         ),
         mesh
     );
-    //mod by alice
-    Info<< "Reading physical velocity field U" << endl;
-    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
-    volVectorField Us
-    (
-        IOobject
-        (
-            "Us",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
 
 //========================
 // drag law modelling
@@ -76,9 +61,8 @@
         mesh
     );
 
-
     //mod by alice
-    Info<< "Reading field phiIB\n" << endl;
+    Info<< "Reading field voidfraction\n" << endl;
     volScalarField voidfraction
     (
         IOobject
@@ -91,6 +75,21 @@
         ),
         mesh
     );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
 //========================
 
 #   include "createPhi.H"
@@ -126,3 +125,5 @@
     );
 
 //===========================
+
+#include "createMRF.H"

applications/solvers/cfdemSolverIBContinuousForcing/Make/files

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

applications/solvers/cfdemSolverIBContinuousForcing/Make/options

@@ -0,0 +1,34 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -ldynamicFvMesh \
+    -ldynamicMesh \
+    -lfvOptions \
+    -lsampling \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
+

applications/solvers/cfdemSolverIBContinuousForcing/UEqn.H

@@ -0,0 +1,19 @@
+fvVectorMatrix UEqn
+(
+    fvm::ddt(voidfractionNext,U) + MRF.DDt(U)
+  + fvm::div(phi, U)
+  + turbulence->divDevReff(U)
+ ==
+    fvOptions(U)
+ +  (lambda*(1-voidfractionNext)/U.mesh().time().deltaT())*(fvc::Sp(1,Us)-fvm::Sp(1,U))
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (piso.momentumPredictor())
+{
+    solve(UEqn == -fvc::grad(p));
+    fvOptions.correct(U);
+}

applications/solvers/cfdemSolverIBContinuousForcing/cfdemSolverIBContinuousForcing.C

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-2015 DCS Computing GmbH,Linz
+                                Copyright (C) 2015-     JKU, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverIBContinuousForcing
+
+Description
+    Transient solver for incompressible flow.
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
+    where additional functionality for CFD-DEM coupling using immersed body
+    (fictitious domain) method and a continuous forcing approach is added.
+Contributions
+    Alice Hager
+    Achuth N. Balachandran Nair
+\*---------------------------------------------------------------------------*/
+
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+
+#include "cfdemCloudIB.H"
+#include "implicitCouple.H"
+
+#include "averagingModel.H"
+#include "regionModel.H"
+#include "voidFractionModel.H"
+
+#include "dynamicFvMesh.H"
+
+#include "cellSet.H"
+
+#include "fvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudIB particleCloud(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        //=== dyM ===================
+        interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
+        mesh.update(); //dyM
+
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        // do particle stuff
+        Info << "- evolve()" << endl;
+        particleCloud.evolve(Us);
+
+        volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
+
+        // Pressure-velocity PISO corrector
+        {
+            MRF.correctBoundaryVelocity(U);
+
+            // Momentum predictor
+            #include "UEqn.H"
+
+            // --- PISO loop
+            while (piso.correct())
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        laminarTransport.correct();
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverIBContinuousForcing/createFields.H

@@ -0,0 +1,143 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading physical velocity field U" << endl;
+    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading the penalization factor field lambda\n" << endl;
+    volScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// drag law modelling
+//========================
+
+    Info<< "\nCreating dummy density field rho = 1\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+
+    Info<< "Reading field phiIB\n" << endl;
+    volScalarField phiIB
+    (
+        IOobject
+        (
+            "phiIB",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    //mod by alice
+    Info<< "Reading field voidfraction\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+
+#   include "createPhi.H"
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+//=== dyM ===================
+
+    Info<< "Reading field interFace\n" << endl;
+    volScalarField interFace
+    (
+        IOobject
+        (
+            "interFace",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        //dimensionedScalar("0", dimensionSet(0, -1, 0, 0, 0), 0.0)
+        dimensionedScalar("0", dimensionSet(0, 0, 0, 0, 0), 0.0)
+    );
+
+//===========================
+
+#include "createMRF.H"

applications/solvers/cfdemSolverIBContinuousForcing/pEqn.H

@@ -0,0 +1,35 @@
+volScalarField rUA = 1.0/UEqn.A();
+surfaceScalarField rUAf(fvc::interpolate(rUA));
+
+U = rUA*UEqn.H();
+
+phi = (fvc::interpolate(U) & mesh.Sf())
+    + rUAf*fvc::ddtCorr(U, phi);            // Is there additional flux term due to the particle presence?
+
+MRF.makeRelative(phi);
+
+adjustPhi(phi, U, p);
+
+while (piso.correctNonOrthogonal())
+{
+    // Pressure corrector
+
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rUA, p) == fvc::div(phi) + particleCloud.ddtVoidfraction()
+    );
+
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+    if (piso.finalNonOrthogonalIter())
+    {
+        phi -= pEqn.flux();
+    }
+}
+
+#include "continuityErrs.H"
+U -= rUA*fvc::grad(p); // should we add a pressure correction?
+U.correctBoundaryConditions();
+fvOptions.correct(U);

applications/solvers/cfdemSolverMultiphaseScalar/Allwclean

@@ -2,7 +2,7 @@
 cd ${0%/*} || exit 1    # Run from this directory
 set -x
 
-wclean libso multiphaseMixture
+wclean libso multiphaseMixtureScalar
 wclean
 
 #------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphaseScalar/Allwmake

@@ -6,7 +6,7 @@ targetType=libso
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 set -x
 
-wmake $targetType multiphaseMixture
+wmake $targetType multiphaseMixtureScalar
 wmake
 
 #------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphaseScalar/Make/options

@@ -6,7 +6,7 @@ include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 EXE_INC = \
      $(PFLAGS) \
     -I$(CFDEM_OFVERSION_DIR) \
-    -ImultiphaseMixture/lnInclude \
+    -ImultiphaseMixtureScalar/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \

applications/solvers/cfdemSolverMultiphaseScalar/cfdemSolverMultiphaseScalar.C

@@ -30,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "multiphaseMixture.H"
+#include "multiphaseMixtureScalar.H"
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"

applications/solvers/cfdemSolverMultiphaseScalar/createFields.H

@@ -88,7 +88,7 @@ surfaceScalarField phi
      linearInterpolate(U*voidfraction) & mesh.Sf()
  );
 
-multiphaseMixture mixture(U, phi, voidfraction);
+multiphaseMixtureScalar mixture(U, phi, voidfraction);
 
 // Need to store rho for ddt(rho, U)
 volScalarField rho

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/Make/files

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

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -26,7 +26,7 @@ Class
 
 Description
     Contact-angle boundary condition for multi-phase interface-capturing
-    simulations.  Used in conjuction with multiphaseMixture.
+    simulations.  Used in conjuction with multiphaseMixtureScalar.
 
 SourceFiles
     alphaContactAngleFvPatchScalarField.C
@@ -37,7 +37,7 @@ SourceFiles
 #define alphaContactAngleFvPatchScalarField_H
 
 #include "zeroGradientFvPatchFields.H"
-#include "multiphaseMixture.H"
+#include "multiphaseMixtureScalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -117,8 +117,8 @@ public:
     typedef HashTable
     <
         interfaceThetaProps,
-        multiphaseMixture::interfacePair,
-        multiphaseMixture::interfacePair::hash
+        multiphaseMixtureScalar::interfacePair,
+        multiphaseMixtureScalar::interfacePair::hash
     > thetaPropsTable;
 
 

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/multiphaseMixtureScalar.C

@@ -18,7 +18,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "multiphaseMixture.H"
+#include "multiphaseMixtureScalar.H"
 #include "alphaContactAngleFvPatchScalarField.H"
 #include "Time.H"
 #include "subCycle.H"
@@ -31,13 +31,13 @@ License
 
 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
 
-const Foam::scalar Foam::multiphaseMixture::convertToRad =
+const Foam::scalar Foam::multiphaseMixtureScalar::convertToRad =
     Foam::constant::mathematical::pi/180.0;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::multiphaseMixture::calcAlphas()
+void Foam::multiphaseMixtureScalar::calcAlphas()
 {
     scalar level = 0.0;
     alphas_ == 0.0;
@@ -51,7 +51,7 @@ void Foam::multiphaseMixture::calcAlphas()
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::calcNu() const
+Foam::multiphaseMixtureScalar::calcNu() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -74,7 +74,7 @@ Foam::multiphaseMixture::calcNu() const
 }
 
 Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseMixture::calcStf() const
+Foam::multiphaseMixtureScalar::calcStf() const
 {
     tmp<surfaceScalarField> tstf
     (
@@ -134,7 +134,7 @@ Foam::multiphaseMixture::calcStf() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::multiphaseMixture::multiphaseMixture
+Foam::multiphaseMixtureScalar::multiphaseMixtureScalar
 (
     const volVectorField& U,
     const surfaceScalarField& phi,
@@ -230,7 +230,7 @@ Foam::multiphaseMixture::multiphaseMixture
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::rho() const
+Foam::multiphaseMixtureScalar::rho() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -247,7 +247,7 @@ Foam::multiphaseMixture::rho() const
 
 
 Foam::tmp<Foam::scalarField>
-Foam::multiphaseMixture::rho(const label patchi) const
+Foam::multiphaseMixtureScalar::rho(const label patchi) const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -264,9 +264,9 @@ Foam::multiphaseMixture::rho(const label patchi) const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::mu() const
+Foam::multiphaseMixtureScalar::mu() const
 {
-    Info << "In multiphasemixture mu()" << endl;
+    Info << "In multiphaseMixtureScalar mu()" << endl;
     return rho()*nu();
 //    PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -283,7 +283,7 @@ Foam::multiphaseMixture::mu() const
 
 
 Foam::tmp<Foam::scalarField>
-Foam::multiphaseMixture::mu(const label patchi) const
+Foam::multiphaseMixtureScalar::mu(const label patchi) const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -306,7 +306,7 @@ Foam::multiphaseMixture::mu(const label patchi) const
 
 
 Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseMixture::muf() const
+Foam::multiphaseMixtureScalar::muf() const
 {
 
     return nuf()*fvc::interpolate(rho());
@@ -327,13 +327,13 @@ Foam::multiphaseMixture::muf() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::nu() const
+Foam::multiphaseMixtureScalar::nu() const
 {
     return nu_;
 }
 
 Foam::tmp<Foam::scalarField>
-Foam::multiphaseMixture::nu(const label patchi) const
+Foam::multiphaseMixtureScalar::nu(const label patchi) const
 {
     //return nu_.boundaryField()[patchi];
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
@@ -355,7 +355,7 @@ Foam::multiphaseMixture::nu(const label patchi) const
 
 
 Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseMixture::nuf() const
+Foam::multiphaseMixtureScalar::nuf() const
 {
     //return muf()/fvc::interpolate(rho());
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
@@ -374,7 +374,7 @@ Foam::multiphaseMixture::nuf() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::Cp() const
+Foam::multiphaseMixtureScalar::Cp() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -395,7 +395,7 @@ Foam::multiphaseMixture::Cp() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::kf() const
+Foam::multiphaseMixtureScalar::kf() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -417,7 +417,7 @@ Foam::multiphaseMixture::kf() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::D() const
+Foam::multiphaseMixtureScalar::D() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -439,7 +439,7 @@ Foam::multiphaseMixture::D() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::Cs() const
+Foam::multiphaseMixtureScalar::Cs() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -456,7 +456,7 @@ Foam::multiphaseMixture::Cs() const
 }
 
 Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseMixture::diffusionCorrection() const
+Foam::multiphaseMixtureScalar::diffusionCorrection() const
 {
 
     surfaceScalarField numerator
@@ -517,7 +517,7 @@ Foam::multiphaseMixture::diffusionCorrection() const
     return correction;
 }
 
-void Foam::multiphaseMixture::solve()
+void Foam::multiphaseMixtureScalar::solve()
 {
     correct();
 
@@ -570,7 +570,7 @@ void Foam::multiphaseMixture::solve()
 }
 
 
-void Foam::multiphaseMixture::correct()
+void Foam::multiphaseMixtureScalar::correct()
 {
     forAllIter(PtrDictionary<phase>, phases_, iter)
     {
@@ -579,7 +579,7 @@ void Foam::multiphaseMixture::correct()
 }
 
 
-Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
+Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureScalar::nHatfv
 (
     const volScalarField& alpha1,
     const volScalarField& alpha2
@@ -605,7 +605,7 @@ Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
 }
 
 
-Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
+Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixtureScalar::nHatf
 (
     const volScalarField& alpha1,
     const volScalarField& alpha2
@@ -622,7 +622,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
 // 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
+void Foam::multiphaseMixtureScalar::correctContactAngle
 (
     const phase& alpha1,
     const phase& alpha2,
@@ -726,7 +726,7 @@ void Foam::multiphaseMixture::correctContactAngle
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureScalar::K
 (
     const phase& alpha1,
     const phase& alpha2
@@ -742,7 +742,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::nearInterface() const
+Foam::multiphaseMixtureScalar::nearInterface() const
 {
     tmp<volScalarField> tnearInt
     (
@@ -768,7 +768,7 @@ Foam::multiphaseMixture::nearInterface() const
 }
 
 
-void Foam::multiphaseMixture::solveAlphas
+void Foam::multiphaseMixtureScalar::solveAlphas
 (
     const scalar cAlpha
 )
@@ -901,7 +901,7 @@ void Foam::multiphaseMixture::solveAlphas
 }
 
 
-bool Foam::multiphaseMixture::read()
+bool Foam::multiphaseMixtureScalar::read()
 {
     if (transportModel::read())
     {

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/multiphaseMixtureScalar.H

@@ -17,10 +17,10 @@ License
     Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
 
 Class
-    multiphaseMixture
+    multiphaseMixtureScalar
 
 Description
-    This class is based on the OpenFOAM(R) Foam::multiphaseMixture class,
+    This class is based on the OpenFOAM(R) Foam::multiphaseMixtureScalar 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
@@ -33,11 +33,11 @@ Description
     between each phase-pair.
 
 SourceFiles
-    multiphaseMixture.C
+    multiphaseMixtureScalar.C
 \*---------------------------------------------------------------------------*/
 
-#ifndef multiphaseMixture_H
-#define multiphaseMixture_H
+#ifndef multiphaseMixtureScalar_H
+#define multiphaseMixtureScalar_H
 
 #include "incompressible/transportModel/transportModel.H"
 #include "IOdictionary.H"
@@ -52,10 +52,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class multiphaseMixture Declaration
+                      Class multiphaseMixtureScalar Declaration
 \*---------------------------------------------------------------------------*/
 
-class multiphaseMixture
+class multiphaseMixtureScalar
 :
     public IOdictionary,
     public transportModel
@@ -191,7 +191,7 @@ public:
     // Constructors
 
         //- Construct from components
-        multiphaseMixture
+        multiphaseMixtureScalar
         (
             const volVectorField& U,
             const surfaceScalarField& phi,
@@ -200,7 +200,7 @@ public:
 
 
     //- Destructor
-    virtual ~multiphaseMixture()
+    virtual ~multiphaseMixtureScalar()
     {}
 
 

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/phase/phase.H

@@ -26,7 +26,7 @@ Class
 
 Description
     Single incompressible phase derived from the phase-fraction.
-    Used as part of the multiPhaseMixture for interface-capturing multi-phase
+    Used as part of the multiphaseMixtureScalar for interface-capturing multi-phase
     simulations.
 
 SourceFiles

applications/solvers/cfdemSolverPisoFreeStreaming/Make/options

@@ -12,7 +12,8 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -51,8 +51,4 @@
     thermo.correct();
 
     Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
 }

applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H

@@ -55,8 +55,4 @@ fvScalarMatrix EEqn
     Info << "Cpv :" << max(Cpv).value() << " " << min(Cpv).value() << endl;
     Info<<  "T max/min : " << max(T).value() << " " << min(T).value() << endl;
     Info << "he max/min : " << max(he).value() << "  " << min(he).value() << endl;
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
 }

applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H

@@ -60,13 +60,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     if (propagateInertSpecie)
     {
-        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + VSMALL);
+        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + ROOTVSMALL);
         forAll(Y,i)
         {
             if (i!=inertIndex)
             {
                 volScalarField& Yi = Y[i];
-                Yi = Yi/(Yt+VSMALL);
+                Yi = Yi/(Yt+ROOTVSMALL);
             }
         }
     }

applications/solvers/cfdemSolverRhoPimpleChem/createFields.H

@@ -303,18 +303,4 @@
         mesh,
         dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0)
     );
-
-    volScalarField dSauter
-    (
-        IOobject
-        (
-            "dSauter",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero",dimensionSet(0, 1, 0, 0, 0,0,0),0)
-    );
 //===============================

applications/solvers/rcfdemSolverCoupledHeattransfer/TEqn.H

@@ -1,4 +1,4 @@
-    volScalarField rhoeps = rhoRec*voidfractionRec;
+    rhoeps = rhoRec*voidfractionRec;
 
     particleCloud.energyContributions(Qsource);
 
@@ -23,10 +23,17 @@
         fvOptions(rhoeps, T)    // no fvOptions support yet
     );
 
-    fvOptions.constrain(TEqn); // no fvOptions support yet
+ //   fvOptions.constrain(TEqn); // no fvOptions support yet
+
+    TEqn.relax();
 
     TEqn.solve();
 
+    T = max(T, TMin);
+    T = min(T, TMax);
+
+    Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
+
     particleCloud.clockM().start(31,"postFlow");
     counter++;
 

applications/solvers/rcfdemSolverCoupledHeattransfer/createFields.H

@@ -73,6 +73,9 @@
         dimensionedVector("zero", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
+    volScalarField rhoeps("rhoeps", rhoRec*voidfractionRec);
+    rhoeps.oldTime(); // switch on saving old time
+
     // heat transfer fields
     Info << "\nCreating heat transfer fields.\n" << endl;
 
@@ -228,3 +231,25 @@
         )
     );
     weightDict.add("weights",scalarList(1,1.0));
+
+    dimensionedScalar TMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "TMax",
+            transportProps,
+            dimTemperature,
+            GREAT
+        )
+    );
+
+    dimensionedScalar TMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "TMin",
+            transportProps,
+            dimTemperature,
+            0.0
+        )
+    );

applications/solvers/rcfdemSolverCoupledHeattransfer/rcfdemSolverCoupledHeattransfer.C

@@ -98,7 +98,7 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().start(26,"Flow");
         #include "updateRho.H"
-        #include "TEqImp.H"
+        #include "TEqn.H"
         particleCloud.clockM().stop("Flow");
 
         stepCounter++;

applications/solvers/rcfdemSolverCoupledHeattransfer/updateRho.H

@@ -1 +1,2 @@
-rhoRec = pRec / (T * R);
+dimensionedScalar Tave = T.weightedAverage(voidfractionRec);
+rhoRec = pRec / (Tave * R);

applications/solvers/rcfdemSolverForcedTracers/Make/options

@@ -12,6 +12,7 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
     -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
@@ -24,4 +25,4 @@ EXE_LIBS = \
     -lfvOptions \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
-$(CFDEM_ADD_LIBS)
+     $(CFDEM_ADD_LIBS)

applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H

@@ -60,9 +60,4 @@
     thermo.correct();
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
-
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
 }

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H

@@ -59,9 +59,4 @@
     thermo.correct();
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
-
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
 }

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options

@@ -27,6 +27,7 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/rctfSpeciesTransport/Make/options

@@ -12,6 +12,7 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
     -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/utilities/displacementField/displacementField.C

@@ -37,11 +37,19 @@ Application
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 void findPairs(labelList &, labelList &, labelPairList &);
 void findPairsUnordered(labelList &, labelList &, labelPairList &);
-void fillEmptyCells(fvMesh &, label, label, labelList &, volVectorField &, volVectorField &, scalarList &, volVectorField &, volVectorField &, bool, scalar);
-void nearestNeighborCells(fvMesh &, label, label, label, labelList &, labelList &);
-void normalizeFields(labelList &, volVectorField &, volVectorField &);
-void readDump(std::string, labelList &, vectorList &);
+void fillEmptyCells(fvMesh &, label, label, scalarList &, volVectorField &, volVectorField &, scalarList &, volVectorField &, volVectorField &, bool, scalar);
+void nearestNeighborCells(fvMesh &, label, label, label, scalarList &, labelList &);
+void normalizeFields(scalarList &, volVectorField &, volVectorField &);
+void readDump(std::string, labelList &, scalarList &, vectorList &);
 scalar weightFun(scalar);
+label maxNumParticles = 1000000;
+scalar minVol = 1e-12;
+scalar Pi43 = 4.1888;
+label posIndex = -1;
+label posRad = -1;
+label posX = -1;
+label posY = -1;
+label posZ = -1;
 
 int main(int argc, char *argv[])
 {
@@ -88,6 +96,11 @@ int main(int argc, char *argv[])
     label dumpIndexDisplacementIncrement(readLabel(displacementProperties.lookup("dumpIndexDisplacementIncrement")));
     label nNeighMin(readLabel(displacementProperties.lookup("nNeighMin")));
     label maxSearchLayers(displacementProperties.lookupOrDefault<label>("maxSearchLayers",0));
+    posIndex = readLabel(displacementProperties.lookup("posIndex"));
+    posRad = readLabel(displacementProperties.lookup("posRad"));
+    posX = readLabel(displacementProperties.lookup("posX"));
+    posY = readLabel(displacementProperties.lookup("posY"));
+    posZ = readLabel(displacementProperties.lookup("posZ"));
     scalar timePerInputStep(readScalar(displacementProperties.lookup("timePerInputStep")));
     scalar timePerDisplacementStep(readScalar(displacementProperties.lookup("timePerDisplacementStep")));
     scalar startTime(readScalar(displacementProperties.lookup("startTime")));
@@ -95,7 +108,7 @@ int main(int argc, char *argv[])
     std::string fileext=string(displacementProperties.lookupOrDefault<string>("fileextension",""));
     bool interpolate=bool(displacementProperties.lookupOrDefault<bool>("fillEmptyCells",true));
     bool averageMode=bool(displacementProperties.lookupOrDefault<bool>("averageMode",false));
-    
+
     volVectorField defaultUs
     (
         IOobject
@@ -110,11 +123,11 @@ int main(int argc, char *argv[])
         dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
     );
 
-    volVectorField defaultUsDirectedVariance
+    volVectorField defaultUsDirectedStdDev
     (
         IOobject
         (
-            "defaultUDispDirectedVariance",
+            "defaultUDispDirectedStdDev",
             runTime.timeName(),
             mesh,
             IOobject::READ_IF_PRESENT,
@@ -172,11 +185,11 @@ int main(int argc, char *argv[])
         dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
     );
 
-    volVectorField UsDirectedVariance
+    volVectorField UsDirectedStdDev
     (
         IOobject
         (
-            "UDispDirectedVariance",
+            "UDispDirectedStdDev",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -186,7 +199,7 @@ int main(int argc, char *argv[])
         dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
     );
 
-    labelList particlesInCell(mesh.nCells(), 0);
+    scalarList particleVolInCell(mesh.nCells(), 0.0);
 
     scalar currTime=startTime + thisProc * timePerInputStep;
     label timeIndex=thisProc;
@@ -196,6 +209,7 @@ int main(int argc, char *argv[])
         runTime.setTime(currTime,timeIndex);
         // read dump files and check which particle indices are present in both
         labelList indices1, indices2;
+        scalarList radii1, radii2;
         vectorList positions1, positions2;
 
         std::stringstream ss;
@@ -209,13 +223,13 @@ int main(int argc, char *argv[])
         {
             if (averageMode)
             {
-                normalizeFields(particlesInCell, Us, UsDirectedVariance);
-                fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
+                normalizeFields(particleVolInCell, Us, UsDirectedStdDev);
+                fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particleVolInCell,Us,UsDirectedStdDev,boundaries,defaultUs,defaultUsDirectedStdDev,interpolate,timePerDisplacementStep);
 
                 Us /= timePerDisplacementStep;
-                UsDirectedVariance /= timePerDisplacementStep;
+                UsDirectedStdDev /= timePerDisplacementStep;
                 Us.write();
-                UsDirectedVariance.write();
+                UsDirectedStdDev.write();
             }
             break;
         }
@@ -225,8 +239,8 @@ int main(int argc, char *argv[])
         Info << "\t" << filename2 << endl;
         Info << "corresponding to time = " << currTime << "." << endl;
 
-        readDump(filename1, indices1, positions1);
-        readDump(filename2, indices2, positions2);
+        readDump(filename1, indices1, radii1, positions1);
+        readDump(filename2, indices2, radii2, positions2);
 
         labelPairList pairs;
         findPairs(indices1,indices2,pairs);
@@ -234,15 +248,16 @@ int main(int argc, char *argv[])
         // average particle displacements and their variance
         Info << "Binning particle displacements on mesh." << endl;
         vector position, displacement;
+        scalar radius, volume;
         label line1, line2;
         label cellI;
 
         if (!averageMode)
         {
             Us *= 0.0;
-            UsDirectedVariance *= 0.0;
-            particlesInCell.clear();
-            particlesInCell.setSize(mesh.nCells(), 0);
+            UsDirectedStdDev *= 0.0;
+            particleVolInCell.clear();
+            particleVolInCell.setSize(mesh.nCells(), 0);
         }
 
         for (label partI = 0; partI < pairs.size(); partI++)
@@ -250,27 +265,29 @@ int main(int argc, char *argv[])
             line1 = pairs[partI].first();
             line2 = pairs[partI].second();
             position = positions1[line1];
-            displacement = positions2[line2] - positions1[line1];
             cellI = mesh.findCell(position);
             if (cellI < 0)  continue;
-            particlesInCell[cellI] += 1;
-            Us[cellI] += displacement;
+            displacement = positions2[line2] - positions1[line1];
+            radius = radii1[line1];
+            volume = Pi43 * radius * radius * radius;
+            particleVolInCell[cellI] += volume;
+            Us[cellI] += displacement*volume;
 
             for (label comp=0;comp<3;comp++)
             {
-                UsDirectedVariance[cellI].component(comp) += displacement.component(comp)*displacement.component(comp);
+                UsDirectedStdDev[cellI].component(comp) += displacement.component(comp)*displacement.component(comp)*volume;
             }
         }
 
         if (!averageMode)
         {
-            normalizeFields(particlesInCell, Us, UsDirectedVariance);
-            fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
+            normalizeFields(particleVolInCell, Us, UsDirectedStdDev);
+            fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particleVolInCell,Us,UsDirectedStdDev,boundaries,defaultUs,defaultUsDirectedStdDev,interpolate,timePerDisplacementStep);
 
             Us /= timePerDisplacementStep;
-            UsDirectedVariance /= timePerDisplacementStep;
+            UsDirectedStdDev /= timePerDisplacementStep;
             Us.write();
-            UsDirectedVariance.write();
+            UsDirectedStdDev.write();
         }
 
         if (averageMode && monitorProbes)
@@ -280,7 +297,7 @@ int main(int argc, char *argv[])
             {
                 vector pos = probePoints[p];
                 label cellP = mesh.findCell(pos);
-                monitoringDataFile << " " << particlesInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedVariance[cellP]/(timePerDisplacementStep*timePerDisplacementStep);
+                monitoringDataFile << " " << particleVolInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedStdDev[cellP]/(timePerDisplacementStep*timePerDisplacementStep);
             }
             monitoringDataFile << endl;
         }
@@ -293,30 +310,69 @@ int main(int argc, char *argv[])
     return 0;
 }
 
-void readDump(std::string filename, labelList &indices, vectorList &positions)
+void readDump(std::string filename, labelList &indices, scalarList &radii, vectorList &positions)
 {
     #include <fstream>
 
     const label leadingLines = 9;
     label lineCounter = 0;
     label partIndex;
-    scalar x, y, z;
+    scalar r = 1.0, x = 0.0, y = 0.0, z = 0.0;
 
     indices.clear();
+    radii.clear();
     positions.clear();
 
+    indices.setSize(maxNumParticles);
+    radii.setSize(maxNumParticles);
+    positions.setSize(maxNumParticles);
+
     std::ifstream file(filename);
     std::string str;
+    std::string word;
+    label wordcounter;
     while (std::getline(file, str))
     {
         if (lineCounter >= leadingLines)
         {
-            sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
-            indices.append(partIndex);
-            positions.append(vector(x,y,z));
+            std::istringstream ss(str);
+            wordcounter = 0;
+            while (ss >> word)
+            {
+                if (wordcounter == posIndex)
+                {
+                    partIndex = stoi(word);
+                }
+                else if (wordcounter == posRad)
+                {
+                    r = stod(word);
+                }
+                else if (wordcounter == posX)
+                {
+                    x = stod(word);
+                }
+                else if (wordcounter == posY)
+                {
+                    y = stod(word);
+                }
+                else if (wordcounter == posZ)
+                {
+                    z = stod(word);
+                }
+                wordcounter++;
+            }
+//            sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
+            indices[lineCounter-leadingLines] = partIndex;
+            radii[lineCounter-leadingLines] = r;
+            positions[lineCounter-leadingLines] = vector(x,y,z);
         }
         lineCounter++;
     }
+
+    label readLines = lineCounter - leadingLines;
+    indices.resize(readLines);
+    radii.resize(readLines);
+    positions.resize(readLines);
 }
 
 void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
@@ -324,6 +380,10 @@ void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
     // remove all entries from first list if they are not present in second list
     // this assumes ordered entries
 
+    pairs.clear();
+    pairs.setSize(maxNumParticles);
+    label pairCounter = 0;
+
     if (indices2.size() == 0) return;
 
     for (label i=0;i<indices1.size();i++)
@@ -339,18 +399,23 @@ void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
             else if (indices2[jmid] < index1) j1 = jmid;
             else
             {
-                pairs.append(labelPair(i,jmid));
+                pairs[pairCounter]=labelPair(i,jmid);
+                pairCounter++;
                 break;
             }
             if (j2-j1 == 1) break;
         }
     }
+    pairs.resize(pairCounter);
     Info << "findPairs: " << pairs.size() << " pairs found." << endl;
 }
 
 void findPairsUnordered(labelList &indices1, labelList &indices2, labelPairList &pairs)
 {
     // remove all entries from first list if they are not present in second list
+    pairs.clear();
+    pairs.setSize(maxNumParticles);
+    label pairCounter = 0;
 
     for (label i=0;i<indices1.size();i++)
     {
@@ -358,15 +423,17 @@ void findPairsUnordered(labelList &indices1, labelList &indices2, labelPairList
         {
             if (indices1[i] == indices2[j])
             {
-                pairs.append(labelPair(i,j));
+                pairs[pairCounter]=labelPair(i,j);
+                pairCounter++;
                 break;
             }
         }
     }
+    pairs.resize(pairCounter);
     Info << "findPairs: " << pairs.size() << " pairs found." << endl;
 }
 
-void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelList &particlesInCell, volVectorField &Us, volVectorField& UsDirectedVariance,scalarList& boundaries, volVectorField &defaultUs, volVectorField &defaultUsDirectedVariance, bool interpolate, scalar dt)
+void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, scalarList &particleVolInCell, volVectorField &Us, volVectorField& UsDirectedStdDev,scalarList& boundaries, volVectorField &defaultUs, volVectorField &defaultUsDirectedStdDev, bool interpolate, scalar dt)
 {
     labelList neighborsWithValues;
     scalar neighborSqrDistance;
@@ -377,7 +444,7 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
     Info << "Filling empty cells." << endl;
     forAll(mesh.C(), cellI)
     {
-        if (particlesInCell[cellI] > 0) continue;
+        if (particleVolInCell[cellI] > minVol) continue;
 
         vector position = mesh.C()[cellI];
         label outsideBox = 0;
@@ -388,11 +455,11 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
         if (outsideBox > 0 || !interpolate)
         {
             Us[cellI] = defaultUs[cellI]*dt;
-            UsDirectedVariance[cellI] = defaultUsDirectedVariance[cellI]*dt;
+            UsDirectedStdDev[cellI] = defaultUsDirectedStdDev[cellI]*dt;
             continue;
         }
 
-        nearestNeighborCells(mesh, cellI, nNeighMin, maxSearchLayers, particlesInCell, neighborsWithValues);
+        nearestNeighborCells(mesh, cellI, nNeighMin, maxSearchLayers, particleVolInCell, neighborsWithValues);
         weightSum = 0.0;
         weights.clear();
         for (label neighI=0; neighI<neighborsWithValues.size(); neighI++)
@@ -406,13 +473,13 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
         {
             weight = weights[neighI]/weightSum;
             Us[cellI] += weight*Us[neighborsWithValues[neighI]];
-            UsDirectedVariance[cellI] += weight*UsDirectedVariance[neighborsWithValues[neighI]];
+            UsDirectedStdDev[cellI] += weight*UsDirectedStdDev[neighborsWithValues[neighI]];
         }
-        
+
         if (neighborsWithValues.size() == 0)
         {
             Us[cellI] = defaultUs[cellI]*dt;
-            UsDirectedVariance[cellI] = defaultUsDirectedVariance[cellI]*dt;
+            UsDirectedStdDev[cellI] = defaultUsDirectedStdDev[cellI]*dt;
         }
 
         // make sure no particles are placed outside of domain
@@ -429,7 +496,7 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
     }
 }
 
-void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label maxSearchLayers, labelList &particlesInCell, labelList &neighborsWithValues)
+void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label maxSearchLayers, scalarList &particleVolInCell, labelList &neighborsWithValues)
 {
     label numSearchLayers = 0;
     std::set<label> neighbors;
@@ -455,11 +522,11 @@ void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label ma
                 {
                     newNeighbors.insert(adj);
                     neighbors.insert(adj);
-                    if (particlesInCell[adj] > 0) neighborsWithValues.append(adj);
+                    if (particleVolInCell[adj] > minVol) neighborsWithValues.append(adj);
                 }
             }
         }
-        
+
         numSearchLayers++;
         if (numSearchLayers > maxSearchLayers && maxSearchLayers > 0) return;
 
@@ -470,18 +537,18 @@ void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label ma
     }
 }
 
-void normalizeFields(labelList& particlesInCell, volVectorField& Us, volVectorField & UsDirectedVariance)
+void normalizeFields(scalarList& particleVolInCell, volVectorField& Us, volVectorField & UsDirectedStdDev)
 {
-    for (label cellJ = 0; cellJ<particlesInCell.size(); cellJ++)
+    for (label cellJ = 0; cellJ<particleVolInCell.size(); cellJ++)
     {
-        if (particlesInCell[cellJ] > 0)
+        if (particleVolInCell[cellJ] > minVol)
         {
-            Us[cellJ] /= particlesInCell[cellJ];
-            UsDirectedVariance[cellJ] /= particlesInCell[cellJ];
+            Us[cellJ] /= particleVolInCell[cellJ];
+            UsDirectedStdDev[cellJ] /= particleVolInCell[cellJ];
             for (label comp=0;comp<3;comp++)
             {
-                UsDirectedVariance[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
-                if (UsDirectedVariance[cellJ].component(comp) > 0) UsDirectedVariance[cellJ].component(comp) = Foam::sqrt(UsDirectedVariance[cellJ].component(comp));
+                UsDirectedStdDev[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
+                if (UsDirectedStdDev[cellJ].component(comp) > 0) UsDirectedStdDev[cellJ].component(comp) = Foam::sqrt(UsDirectedStdDev[cellJ].component(comp));
             }
         }
     }

doc/CFDEMcoupling_models.txt

@@ -124,6 +124,7 @@ particleDeformation,
 potentialRelaxation,
 "surfaceTensionForce"_forceModel_surfaceTensionForce.html,
 terminalVelocity,
+"transferFluidProperties"_forceModel_transferFluidProperties.html,
 turbulentDispersion,
 turbulentVelocityFluctuations,
 "virtualMassForce"_forceModel_virtualMassForce.html,

doc/CFDEMcoupling_solvers.txt

@@ -10,6 +10,7 @@
 This section lists all CFDEMcoupling solvers alphabetically.
 
 "cfdemSolverIB"_cfdemSolverIB.html,
+"cfdemSolverIBContinuousForcing"_cfdemSolverIBContinuousForcing.html,
 "cfdemSolverMultiphase"_cfdemSolverMultiphase.html,
 "cfdemSolverMultiphaseScalar"_cfdemSolverMultiphaseScalar.html,
 "cfdemSolverPiso"_cfdemSolverPiso.html,

doc/IOModel_basicIO.txt

@@ -17,7 +17,8 @@ IOModel basicIO; :pre
 [Examples:]
 
 IOModel basicIO;
-serialOutput; :pre
+serialOutput;
+cartesianOutput; :pre
 
 [Description:]
 
@@ -29,6 +30,11 @@ Using the keyword {serialOutput;} in the
 the IO is serial to the directory {$casePath/CFD/lagrangian}. In this case
 only the data on processor 0 is written!
 
+Using the keyword {cartesianOutput;} in the
+"couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties dictionary,
+particle positions will be written as cartesian coordinates instead of
+barycentric coordinates.
+
 Data is written every write time of the CFD simulation.
 
 [Restrictions:]

doc/cfdemSolverIBContinuousForcing.txt

@@ -0,0 +1,82 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+cfdemSolverIBContinuousForcing command :h3
+
+[Description:]
+
+"cfdemSolverIBContinuousForcing" is a coupled CFD-DEM solver using CFDEMcoupling,
+an open-source parallel coupled CFD-DEM framework, for calculating the dynamics
+between immersed bodies and the surrounding fluid. Being an implementation of a
+continuous forcing immersed boundary method it allows tackling problems where
+the body diameter exceeds the maximal size of a fluid cell.
+
+<!-- HTML_ONLY -->
+Using the toolbox of OpenFOAM&reg;(*) the governing equations of the fluid are
+computed and the corrections of velocity and pressure field with respect to the
+body-movement information, gained from LIGGGHTS, are incorporated.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+Using the toolbox of OpenFOAM\ |reg|\ (*) the governing equations of the fluid are
+computed and the corrections of velocity and pressure field with respect to the
+body-movement information, gained from LIGGGHTS, are incorporated.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+
+The code of this solver was contributed by A.N. Balachandran Nair, JKU. For more
+details, see "Balachandran Nair et al. (2021)"_#BalachandranNair2021
+
+[Use:]
+
+The solver is realized within the open-source framework CFDEMcoupling. Just as
+for the unresolved CFD-DEM solver cfdemSolverPiso the file
+CFD/constant/couplingProperties contains information about the settings for the
+different models. While IOmodel, DataExchangeModel etc. are applicable for all
+CFDEMcoupling-solvers, special locate-, force- and void fraction models were
+designed for this solver:
+
+"engineSearchIB"_locateModel_engineSearchIB.html,
+"ArchimedesIB"_forceModel_ArchimedesIB.html,
+"ShirgaonkarIB"_forceModel_ShirgaonkarIB.html,
+"IBVoidfraction"_voidFractionModel_IBVoidFraction.html
+
+:line
+
+:link(BalachandranNair2021)
+[(Balachandran Nair, 2021)] Balachandran Nair, A.N., Pirker, S. and Saeedipour, M.,
+"Resolved CFD-DEM simulation of blood flow with a reduced-order RBC model",
+Comp. Part. Mech. (2021)
+
+: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/forceModel_ParmarBassetForce.txt

@@ -23,17 +23,19 @@ ParmarBassetForceProps
     nIntegral           scalar1;
     discretisationOrder scalar2;
     treatForceExplicit  switch1;
-    interpolation       switch2;
+    verbose             switch2;
+    interpolation       switch3;
     smoothingModel      "smoothingModel";
 \} :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
 {Us} = name of the finite volume cell averaged particle velocity field :l
-{scalar1} = number of timesteps considered in the near history :l
-{scalar2} = (1 or 2) discretisation order of the far history differential equations :l
+{scalar1} = number of timesteps considered in the near history (integer > 1):l
+{scalar2} = (optional, default 1) discretisation order of the far history differential equations (1 or 2) :l
 {switch1} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
 {switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
-{smoothingModel} = name of smoothing model for the dU/dt field :l
+{switch3} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{smoothingModel} = name of smoothing model for the dUrel/dt field :l
 :ule
 
 [Examples:]

doc/forceModel_ShirgaonkarIB.txt

@@ -21,12 +21,14 @@ ShirgaonkarIBProps
     velFieldName      "U";
     pressureFieldName "pressure";
     twoDimensional;
+    useTorque;
     treatForceExplicit switch1;
 \} :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
 {pressure} = name of the finite volume pressure field :l
 {twoDimensional} = optional keyword for conducting a two dimensional calculation :l
+{useTorque} =  optional keyword for calculating particle torque :l
 {switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
 :ule
 

doc/forceModel_transferFluidProperties.txt

@@ -0,0 +1,42 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+forceModel transferFluidProperties command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+forceModels
+(
+    transferFluidProperties
+);
+transferFluidPropertiesProps
+\{
+    verbose       switch1;
+    interpolation switch2;
+\} :pre
+
+{switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :ulb,l
+{switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+:ule
+
+[Description:]
+
+This "force model" does not influence the particles or the flow - it transfer to fluid density and (dynamic)
+viscosity from OpenFOAM to LIGGGHTS.
+
+
+[Restrictions:]
+
+This model requires {fix cfd/coupling/fluidproperties} to work.
+
+[Related commands:]
+
+"forceModel"_forceModel.html
+

doc/smoothingModel_constDiffSmoothing.txt

@@ -15,17 +15,22 @@ dictionary.
 smoothingModel constDiffSmoothing;
 constDiffSmoothingProps
 \{
-    lowerLimit                    number1;
-    upperLimit                    number2;
-    smoothingLength               lengthScale;
-    smoothingLengthReferenceField lengthScaleRefField;
+    lowerLimit                        number1;
+    upperLimit                        number2;
+    smoothingLength                   lengthScale;
+    smoothingLengthReference          lengthScaleRef;
+    smoothingLengthFieldName          fieldName1;
+    smoothingLengthReferenceFieldName fieldName2;
     verbose;
 \} :pre
 
 {number1} = scalar fields will be bound to this lower value :ulb,l
 {number2} = scalar fields will be bound to this upper value :l
 {lengthScale} = length scale over which the exchange fields will be smoothed out :l
-{lengthScaleRefField} = length scale over which reference fields (e.g., the average particle velocity) will be smoothed out. Should be always larger than lengthScale. If not specified, will be equal to lengthScale. :l
+{lengthScaleRef} = (optional) length scale over which reference fields (e.g., the average particle velocity) will be smoothed out. Should be always larger than lengthScale. If not specified, will be equal to lengthScale. :l
+{fieldName1} = (optional) name of scalar field to be used as local smoothing length. :l
+{fieldName2} = (optional) name of scalar field to be used as local smoothing length for reference fields. :l
+
 {verbose} = (optional, default false) flag for debugging output :l
 :ule
 
@@ -51,6 +56,11 @@ which these reference fields are not specified. Values calculated in the cells
 e.g. the average particle velocity, which are not specified in all cells in case
 the flow is rather dilute.
 
+Alternative to {smoothingLength} and {smoothingLengthReference}, 
+{smoothingLengthFieldName} and/or {smoothingLengthReferenceFieldName} can be used
+to define spatial variation of the smoothing lengths. Either the scalar or field 
+options must be used, giving both will result in errors.
+
 [Restrictions:]
 
 This model is tested in a limited number of flow situations.

etc/bashrc

@@ -17,7 +17,7 @@
 #------------------------------------------------------------------------------
 
 export CFDEM_PROJECT=CFDEM
-export CFDEM_VERSION=21.11
+export CFDEM_VERSION=23.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 21.11
+setenv CFDEM_VERSION 23.02
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade

etc/functions.sh

@@ -359,7 +359,7 @@ cleanCFDEMcase()
     cd $casePath/CFD
     cleanCase
     rm -r $casePath/DEM/post/*
-    echo "dummyfile" >> $casePath/DEM/post/dummy
+    touch $casePath/DEM/post/.gitignore
     cd $casePath
     echo "done"
 }

etc/library-list.txt

@@ -3,4 +3,4 @@ lagrangian/cfdemParticleComp/dir
 recurrence/dir
 finiteVolume/dir
 ../applications/solvers/cfdemSolverMultiphase/multiphaseMixture/dir
-../applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/dir
+../applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/dir

etc/solver-list.txt

@@ -10,6 +10,7 @@ cfdemSolverPiso/dir
 cfdemSolverRhoPimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
+cfdemSolverIBContinuousForcing/dir
 cfdemSolverRhoPimpleChem/dir
 cfdemSolverMultiphase/dir
 cfdemSolverMultiphaseScalar/dir

etc/tutorial-list.txt

@@ -19,3 +19,5 @@ cfdemSolverIB/twoSpheresGlowinskiMPI/dir
 cfdemSolverPisoScalar/packedBedTemp/dir
 
 cfdemSolverPiso/ErgunTestCG/dir
+
+cfdemSolverIB/falling_sphere_two_way_coupling/dir

src/lagrangian/cfdemParticle/Make/files

@@ -96,6 +96,7 @@ $(forceModels)/potentialRelaxation/potentialRelaxation.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
 $(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
+$(forceModels)/transferFluidProperties/transferFluidProperties.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C
 $(forceModels)/Fines/FanningDynFines.C

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -34,8 +34,8 @@ Description
 #ifndef versionInfo_H
 #define versionInfo_H
 
-word CFDEMversion="PFM 21.11";
-word compatibleLIGGGHTSversion="PFM 21.11";
+word CFDEMversion="PFM 23.02";
+word compatibleLIGGGHTSversion="PFM 23.02";
 word OFversion="6";
 
 Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;

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

@@ -274,6 +274,8 @@ bool cfdemCloudEnergy::evolve
 
 void cfdemCloudEnergy::postFlow()
 {
+    if (ignore()) return;
+
     cfdemCloud::postFlow();
     forAll(energyModel_, modeli)
     {
@@ -285,18 +287,6 @@ void cfdemCloudEnergy::postFlow()
     }
 }
 
-void cfdemCloudEnergy::solve()
-{
-    forAll(energyModel_, modeli)
-    {
-        energyModel_[modeli].solve();
-    }
-    forAll(massTransferModel_, modeli)
-    {
-        massTransferModel_[modeli].solve();
-    }
-}
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -42,7 +42,7 @@ namespace Foam
 class energyModel;
 class massTransferModel;
 class thermCondModel;
-class diffCoeffModel;  
+class diffCoeffModel;
 class chemistryModel;
 /*---------------------------------------------------------------------------*\
                            Class cfdemCloudEnergy Declaration
@@ -142,8 +142,6 @@ public:
 
        void postFlow();
 
-       void solve();
-
 };
 
 

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

@@ -36,7 +36,6 @@ Description
 #include "locateModel.H"
 #include "dataExchangeModel.H"
 #include "IOModel.H"
-#include <mpi.h>
 #include "IOmanip.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,6 +57,7 @@ cfdemCloudIB::cfdemCloudIB
     angularVelocities_(NULL),
     pRefCell_(readLabel(mesh.solutionDict().subDict("PISO").lookup("pRefCell"))),
     pRefValue_(readScalar(mesh.solutionDict().subDict("PISO").lookup("pRefValue"))),
+    DEMTorques_(NULL),
     haveEvolvedOnce_(false),
     skipLagrangeToEulerMapping_(false)
 {
@@ -75,6 +75,7 @@ cfdemCloudIB::cfdemCloudIB
 cfdemCloudIB::~cfdemCloudIB()
 {
     dataExchangeM().destroy(angularVelocities_,3);
+    dataExchangeM().destroy(DEMTorques_,3);
 }
 
 
@@ -82,8 +83,7 @@ cfdemCloudIB::~cfdemCloudIB()
 void cfdemCloudIB::getDEMdata()
 {
     cfdemCloud::getDEMdata();
-    Info << "=== cfdemCloudIB::getDEMdata() === particle rotation not considered in CFD" << endl;
-    //dataExchangeM().getData("omega","vector-atom",angularVelocities_);
+    dataExchangeM().getData("omega","vector-atom",angularVelocities_);
 }
 
 bool cfdemCloudIB::reAllocArrays()
@@ -92,12 +92,24 @@ bool cfdemCloudIB::reAllocArrays()
     {
         // get arrays of new length
         dataExchangeM().allocateArray(angularVelocities_,0.,3);
+        dataExchangeM().allocateArray(DEMTorques_,0.,3);
         return true;
     }
     return false;
 }
 
-bool cfdemCloudIB::evolve()
+void cfdemCloudIB::giveDEMdata()
+{
+    cfdemCloud::giveDEMdata();
+    dataExchangeM().giveData("hdtorque","vector-atom",DEMTorques_);
+}
+
+inline double ** cfdemCloudIB::DEMTorques() const
+{
+    return DEMTorques_;
+}
+
+bool cfdemCloudIB::evolve(volVectorField& Us)
 {
     numberOfParticlesChanged_ = false;
     arraysReallocated_=false;
@@ -140,6 +152,10 @@ bool cfdemCloudIB::evolve()
         for (int i=0;i<nrForceModels();i++) forceM(i).setForce();
         if(verbose_) Info << "setForce done." << endl;
 
+        // set particle velocity field
+        if(verbose_) Info << "- setVelocity(velocities_)" << endl;
+        calcForcingTerm(Us); // only needed for continuous forcing
+
         // write DEM data
         if(verbose_) Info << " -giveDEMdata()" << endl;
         giveDEMdata();
@@ -148,7 +164,8 @@ bool cfdemCloudIB::evolve()
 
         haveEvolvedOnce_=true;
     }
-    Info << "evolve done." << endl;
+
+    if(verbose_) Info << "evolve done." << endl;
 
     //if(verbose_)    #include "debugInfo.H";
 
@@ -158,6 +175,34 @@ bool cfdemCloudIB::evolve()
     return doCouple;
 }
 
+void cfdemCloudIB::calcForcingTerm(volVectorField& Us)
+{
+    label cell = 0;
+    vector uP(0,0,0);
+    vector rVec(0,0,0);
+    vector velRot(0,0,0);
+    vector angVel(0,0,0);
+    Us.primitiveFieldRef() = vector::zero;
+    for (int index = 0; index < numberOfParticles(); ++index)
+    {
+        for (int subCell = 0; subCell < voidFractionM().cellsPerParticle()[index][0]; subCell++)
+        {
+            cell = cellIDs()[index][subCell];
+
+            if (cell >=0)
+            {
+                // calc particle velocity
+                for (int i=0;i<3;i++) rVec[i]=Us.mesh().C()[cell][i]-position(index)[i];
+                for (int i=0;i<3;i++) angVel[i]=angularVelocities()[index][i];
+                velRot=angVel^rVec;
+                for (int i=0;i<3;i++) uP[i] = velocities()[index][i]+velRot[i];
+                Us[cell] = (1-voidfractions_[index][subCell])*uP;
+            }
+        }
+    }
+    Us.correctBoundaryConditions();
+}
+
 void cfdemCloudIB::calcVelocityCorrection
 (
     volScalarField& p,
@@ -190,7 +235,7 @@ void cfdemCloudIB::calcVelocityCorrection
                     for(int i=0;i<3;i++) uParticle[i] = velocities()[index][i]+velRot[i];
 
                     // impose field velocity
-                    U[cellI]=(1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI];
+                    U[cellI]= (1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI];
                 }
             }
         //}
@@ -213,7 +258,7 @@ void cfdemCloudIB::calcVelocityCorrection
     U.correctBoundaryConditions();
 
     // correct the pressure as well
-    p=p+phiIB/U.mesh().time().deltaT();  // do we have to  account for rho here?
+    p=p+phiIB/U.mesh().time().deltaT();  // no need to account for rho since p=(p/rho) in OF
     p.correctBoundaryConditions();
 
     if (couplingProperties_.found("checkinterface"))

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

@@ -62,6 +62,8 @@ protected:
     label pRefCell_;
     scalar pRefValue_;
 
+    double **DEMTorques_;
+
     bool haveEvolvedOnce_;
     bool skipLagrangeToEulerMapping_;
 
@@ -85,7 +87,13 @@ public:
 
         bool reAllocArrays();
 
-        bool evolve();
+        void giveDEMdata();
+
+        inline double ** DEMTorques() const;
+
+        bool evolve(volVectorField&);
+
+        void calcForcingTerm(volVectorField&);
 
         void calcVelocityCorrection(volScalarField&,volVectorField&,volScalarField&,volScalarField&); // this could be moved to an IB mom couple model
 
@@ -96,7 +104,6 @@ public:
         {
             return angularVelocities_;
         }
-
 };
 
 

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

@@ -125,6 +125,14 @@ void IOModel::streamDataToPath(const fileName& path, const double* const* array,
         else fileStream << "{}" << "\n";
         return;
     }
+#if OPENFOAM_VERSION_MAJOR > 4
+    else if (barycentricOutput_ && type == "position")
+    {
+        // Force construction of face-diagonal decomposition before construction
+        // of particle which uses parallel transfers.
+        (void)particleCloud_.mesh().tetBasePtIs();
+    }
+#endif
 
     fileStream << token::BEGIN_LIST << nl;
 
@@ -139,19 +147,24 @@ void IOModel::streamDataToPath(const fileName& path, const double* const* array,
             }
             else if (type == "position")
             {
-#if OPENFOAM_VERSION_MAJOR < 5
-                fileStream << "( " << array[index][0] << " " << array[index][1] << " " << array[index][2] << " ) " << cellIDs[index][0] << nl;
-#else
-                particle part
-                (
-                    particleCloud_.mesh(),
-                    vector(array[index][0],array[index][1],array[index][2]),
-                    cellIDs[index][0]
-                );
+#if OPENFOAM_VERSION_MAJOR > 4
+                if (barycentricOutput_)
+                {
+                    particle part
+                    (
+                        particleCloud_.mesh(),
+                        vector(array[index][0],array[index][1],array[index][2]),
+                        cellIDs[index][0]
+                    );
 
-                part.writePosition(fileStream);
-                fileStream << nl;
+                    part.writePosition(fileStream);
+                    fileStream << nl;
+                }
+                else
 #endif
+                {
+                    fileStream << "( " << array[index][0] << " " << array[index][1] << " " << array[index][2] << " ) " << cellIDs[index][0] << nl;
+                }
             }
             else if (type == "vector")
             {
@@ -179,7 +192,8 @@ IOModel::IOModel
     dict_(dict),
     particleCloud_(sm),
     time_(sm.mesh().time()),
-    parOutput_(true)
+    parOutput_(true),
+    barycentricOutput_(true)
 {
     if (
             particleCloud_.dataExchangeM().type()=="oneWayVTK" ||
@@ -189,6 +203,11 @@ IOModel::IOModel
         parOutput_ = false;
         Warning << "IO model is in serial write mode, only data on proc 0 is written" << endl;
     }
+
+    if (dict_.found("cartesianOutput"))
+    {
+        barycentricOutput_ = false;
+    }
 }
 
 

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

@@ -65,6 +65,8 @@ protected:
 
         bool parOutput_;
 
+        bool barycentricOutput_;
+
 public:
 
     //- Runtime type information

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

@@ -246,8 +246,6 @@ void diffusionCoefficient::execute()
         double**& diffusionCoefficients_ = particleCloud_.getParticlePropertyRef<double**>(diffusantGasNames_[j]);
         particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_);
     }
-
-    Info << "give data done" << endl;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -57,8 +57,10 @@ initMultiLayers::initMultiLayers
     mesh_(sm.mesh()),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
     maxNumLayers_(0),
+    maxNumParticlesPerType_(propsDict_.lookupOrDefault<label>("maxNumParticlesPerType",1000000)),
     numLayers_(propsDict_.lookupOrDefault<labelList>("numLayers",labelList(1,-1))),
     partTypes_(propsDict_.lookupOrDefault<labelList>("partTypes",labelList(1,-1))),
+    searchLayers_(propsDict_.lookupOrDefault<labelList>("searchLayers",labelList(1,-1))),
     relRadiiRegName_(typeName + "relRadii"),
     filepath_(string(propsDict_.lookup("filepath"))),
     initialized_(false)
@@ -107,7 +109,7 @@ bool initMultiLayers::init()
         }
         Info << "\nReading" << endl;
         Info << "\t" << filename << endl;
-        readDump(filename, type, indices, positions, relradii);
+        label numLines = readDump(filename, type, indices, positions, relradii);
 
         Info << "Binning particle displacements on mesh for type " << type << endl;
 
@@ -144,7 +146,7 @@ bool initMultiLayers::init()
         );
 
         label cellI = -1;
-        for (label lineI = 0; lineI < indices.size(); lineI++)
+        for (label lineI = 0; lineI < numLines; lineI++)
         {
             position = positions[lineI];
             relrad = relradii[lineI];
@@ -178,8 +180,8 @@ bool initMultiLayers::init()
             if(cellK >= 0 && partType == type)
             {
                 // look for the nearest occupied cell
-                cellKoccupied = findNearestCellWithValue(cellK, particlesInCell);
                 listIndex = getListIndex(partType);
+                cellKoccupied = findNearestCellWithValue(cellK, particlesInCell,searchLayers_[listIndex]);
                 if (cellKoccupied >= 0)
                 {
                     relRadiiLoopVar = relRadiiField[cellKoccupied];
@@ -215,7 +217,7 @@ void initMultiLayers::execute()
     }
 }
 
-label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &particlesInCell) const
+label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &particlesInCell, label searchLayers) const
 {
     if (particlesInCell[refCell] > 0.5) return refCell;
 
@@ -226,6 +228,8 @@ label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &p
     neighbors.insert(refCell);
     recentNeighbors.insert(refCell);
 
+    label layersSearched = 0;
+
     while(true)
     {
         for (std::set<label>::iterator it=recentNeighbors.begin(); it!=recentNeighbors.end(); ++it)
@@ -250,6 +254,8 @@ label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &p
         recentNeighbors.clear();
         recentNeighbors = newNeighbors;
         newNeighbors.clear();
+        layersSearched++;
+        if (layersSearched >= searchLayers && searchLayers > 0) return -1;
     }
 }
 
@@ -264,7 +270,7 @@ label initMultiLayers::getListIndex(label testElement) const
     return -1;
 }
 
-void initMultiLayers::readDump(std::string filename, label type, labelList &indices, vectorList &positions, vectorList &relradii)
+label initMultiLayers::readDump(std::string filename, label type, labelList &indices, vectorList &positions, vectorList &relradii)
 {
     #include <fstream>
 
@@ -281,6 +287,10 @@ void initMultiLayers::readDump(std::string filename, label type, labelList &indi
     positions.clear();
     relradii.clear();
 
+    indices.setSize(maxNumParticlesPerType_);
+    positions.setSize(maxNumParticlesPerType_);
+    relradii.setSize(maxNumParticlesPerType_);
+
     std::ifstream file(filename);
     std::string str;
     while (std::getline(file, str))
@@ -292,12 +302,17 @@ void initMultiLayers::readDump(std::string filename, label type, labelList &indi
             {
                 FatalError<< "Particle of type " << partType << " detected in " << filename << abort(FatalError);
             }
-            indices.append(partIndex);
-            positions.append(vector(x,y,z));
-            relradii.append(vector(r1,r2,r3));
+
+            indices[lineCounter-leadingLines] = partIndex;
+            positions[lineCounter-leadingLines] = vector(x,y,z);
+            relradii[lineCounter-leadingLines] = vector(r1,r2,r3);
         }
         lineCounter++;
+        if (lineCounter == maxNumParticlesPerType_) break;
     }
+
+    label readLines = lineCounter - leadingLines;
+    return readLines;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -59,11 +59,15 @@ private:
     bool verbose_;
 
     label maxNumLayers_;
+    
+    label maxNumParticlesPerType_;
 
     labelList numLayers_;
 
     labelList partTypes_;
 
+    labelList searchLayers_;
+
     vector defaultRelRadii_;
 
     const word relRadiiRegName_;
@@ -74,9 +78,9 @@ private:
 
     bool init();
 
-    void readDump(std::string, label, labelList &, vectorList &, vectorList &);
+    label readDump(std::string, label, labelList &, vectorList &, vectorList &);
 
-    label findNearestCellWithValue(label, volScalarField &) const;
+    label findNearestCellWithValue(label, volScalarField &, label) const;
 
     label getListIndex(label) const;
 

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

@@ -174,8 +174,6 @@ void massTransferCoeff::execute()
     // give DEM data
     particleCloud_.dataExchangeM().giveData(partNuName_, "scalar-atom", partNu_);
     particleCloud_.dataExchangeM().giveData(partReName_, "scalar-atom", partRe_);
-
-    Info << "give data done" << endl;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -169,8 +169,6 @@ void reactantPerParticle::execute()
 
     // give DEM data
     particleCloud_.dataExchangeM().giveData(partReactantName_, "scalar-atom", reactantPerParticle_);
-
-    Info << "give data done" << endl;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -100,6 +100,7 @@ void One2One::setup
     comm_
   );
 
+  // count number of receives
   ncollected_ = 0;
   int nrequests = 0;
   for (int i = 0; i < nsrc_procs_; i++)
@@ -115,6 +116,18 @@ void One2One::setup
     }
   }
 
+  // count number of sends
+  if (nlocal_ > 0)
+  {
+    for (int i = 0; i < ndst_procs_; i++)
+    {
+      if (dst_procs_[i] != me_)
+      {
+        nrequests++;
+      }
+    }
+  }
+
   if (nrequests > 0)
   {
     request_ = new MPI_Request[nrequests];
@@ -143,23 +156,27 @@ void One2One::exchange(T *&src, T *&dst, int data_length)
     {
       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
+        int tag = (src_procs_[i] << 16) | me_;
+
+        #ifdef O2ODEBUG
+        std::cout<< "[" << me_ << "]"
+                << " RCV " << (i+1)
+                << " of "  << nsrc_procs_
+                << " from: " << src_procs_[i]
+                << " natoms_[src_procs_[i]] " << natoms_[src_procs_[i]]
+                << " data_length " << data_length
+                << " offset " << offset
+                << " tag " << tag
+                << std::endl;
+        #endif
+
         MPI_Irecv
         (
           &dst[offset],
           natoms_[src_procs_[i]]*data_length,
           wrap.mpi_type,
           src_procs_[i],
-          MPI_ANY_TAG,
+          tag,
           comm_,
           &request_[requesti]
         );
@@ -168,16 +185,34 @@ void One2One::exchange(T *&src, T *&dst, int data_length)
       else // data is available on-proc
       {
         offset_local = offset;
+
+        #ifdef O2ODEBUG
+        std::cout<< "[" << me_ << "]"
+                << " RCV " << (i+1)
+                << " of "  << nsrc_procs_
+                << " from: " << src_procs_[i]
+                << " skipped (self)"
+                << " offset_local " << offset_local
+                << std::endl;
+        #endif
       }
     }
+    #ifdef O2ODEBUG
+    else
+    {
+      std::cout<< "[" << me_ << "]"
+              << " RCV " << (i+1)
+              << " of "  << nsrc_procs_
+              << " from: " << src_procs_[i]
+              << " skipped (no data)" 
+              << std::endl;      
+    }
+    #endif
+
     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?
+  // nonblocking sends
   // only do sends if I have particles
   if (nlocal_ > 0)
   {
@@ -185,29 +220,62 @@ void One2One::exchange(T *&src, T *&dst, int data_length)
     {
       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
+        int tag = (me_ << 16) | dst_procs_[i];
+
+        #ifdef O2ODEBUG
+        std::cout<< "[" << me_ << "]"
+                << " SEND to: " << dst_procs_[i]
+                << " nlocal_ " << nlocal_
+                << " data_length " << data_length
+                << " tag " << tag
+                << std::endl;
+        #endif
+
+        MPI_Isend
         (
           src,
           nlocal_*data_length,
           wrap.mpi_type,
           dst_procs_[i],
-          0,
-          comm_
+          tag,
+          comm_,
+          &request_[requesti]
         );
+        requesti++;
       }
+      #ifdef O2ODEBUG
+      else
+      {
+        std::cout<< "[" << me_ << "]"
+                << " SEND to: " << dst_procs_[i]
+                << " skipped (self)" << std::endl;
+      }
+      #endif
     }
   }
+  #ifdef O2ODEBUG
+  else
+  {
+    std::cout<< "[" << me_ << "]" << " SEND skipped (no data)" << std::endl;
+  }
+  #endif
 
   // only wait if requests were actually posted
   if (requesti > 0)
+  {
+    #ifdef O2ODEBUG
+    std::cout<< "[" << me_ << "]"
+            << " WAIT for: " << requesti << " requests" << std::endl;
+    #endif
+
     MPI_Waitall(requesti, request_, status_);
+  }
+  #ifdef O2ODEBUG
+  else
+  {
+    std::cout<< "[" << me_ << "]" << " WAIT skipped (no requests)" << std::endl;
+  }
+  #endif
 
   // copy on-proc data
   if (offset_local > -1)
@@ -223,6 +291,13 @@ void One2One::exchange(T *&src, T *&dst, int data_length)
       dst[locali+offset_local] = src[locali];
     }
   }
+
+  // wait for all procs to complete data exchange
+  #ifdef O2ODEBUG
+  std::cout<< "[" << me_ << "] BARRIER" << std::endl;
+  #endif
+
+  MPI_Barrier(comm_);
 }
 
 template void One2One::exchange<int>(int*&, int*&, int);

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

@@ -183,6 +183,7 @@ void twoWayOne2One::createProcMap()
         point(ligbb[0][0], ligbb[0][1], ligbb[0][2]),
         point(ligbb[1][0], ligbb[1][1], ligbb[1][2])
     );
+    delete [] ligbb;
     if (boundBoxMargin_ > 0.)
     {
         thisLigBox.max() = thisLigBox.max() + boundBoxMargin_ * vector::one;
@@ -230,6 +231,7 @@ twoWayOne2One::~twoWayOne2One()
 {
     delete foam2lig_;
     delete lig2foam_;
+    delete [] lig2foam_mask_;
 
     destroy(lig2foam_ids_);
     destroy(foam2lig_ids_);
@@ -777,10 +779,7 @@ void twoWayOne2One::locateParticles()
         lig2foam_->exchange(my_prev_cell_ids, prev_cell_ids_);
     }
 
-    if (lig2foam_mask_)
-    {
-        delete [] lig2foam_mask_;
-    }
+    delete [] lig2foam_mask_;
     lig2foam_mask_ = new bool[lig2foam_->ncollected_];
 
     DynamicList<label> cellIds;

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

@@ -105,8 +105,6 @@ public:
 
         virtual void postFlow() {}
 
-        virtual void solve() {}
-
         const volScalarField& kf0Field() const;
 
         const volScalarField& CpField() const;
@@ -116,7 +114,6 @@ public:
             return 0.0;
         }
 
-
 };
 
 

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

@@ -283,8 +283,6 @@ void heatTransferInterGrain::calcEnergyContribution()
         if(cellI >= 0)
         {
             voidfraction = voidfraction_[cellI];
-            if (voidfraction < 0.01)
-                voidfraction = 0.01;
 
             partVolume = particleCloud_.particleVolume(index);
             QPartPart = QPartPart_[cellI];
@@ -426,8 +424,6 @@ void heatTransferInterGrain::calcPartThermRad()
 
 scalar heatTransferInterGrain::FE(scalar voidfraction, scalar emissivity, scalar L)
 {
-    if (voidfraction < 1 - SMALL) return 1.0; // limit of voidfraction -> 1.0
-
     scalar B = 1.25 * pow((1-voidfraction)/voidfraction,1.111);
     scalar sqrt_aP = sqrt(1-voidfraction);
     scalar x = 2.0/emissivity - 1.0;

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

@@ -104,7 +104,7 @@ protected:
 
     const scalar TMin = 10.0;
 
-    const scalar voidfracMax = 0.95;
+    const scalar voidfracMax = 0.8; // maximum voidfraction for which B-B correlation is valid; TODO: check value!
 
     const scalar voidfracMin = 0.05;
 

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

@@ -77,7 +77,6 @@ BeetstraDrag::BeetstraDrag
     particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must  be the force
     particleCloud_.probeM().vectorFields_.append("Urel");
     particleCloud_.probeM().scalarFields_.append("Rep");
-    particleCloud_.probeM().scalarFields_.append("betaP");
     particleCloud_.probeM().scalarFields_.append("voidfraction");
     particleCloud_.probeM().writeHeader();
 

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

@@ -34,7 +34,6 @@ Description
 #include "ParmarBassetForce.H"
 #include "addToRunTimeSelectionTable.H"
 #include "smoothingModel.H"
-#include "constDiffSmoothing.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -71,17 +70,18 @@ ParmarBassetForce::ParmarBassetForce
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     UsFieldName_(propsDict_.lookup("granVelFieldName")),
     Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
-    nInt_(readLabel(propsDict_.lookup("nIntegral"))),
-    discOrder_(readLabel(propsDict_.lookup("discretisationOrder"))),
-    nHist_(nInt_+discOrder_+1),
-    FHistSize_(2*discOrder_+1),
-    ddtUrelHist_(nHist_,NULL),                     // UrelHist_[ndt in past][particle ID][dim]
-    rHist_(nHist_,NULL),                           // rHist_[ndt in past][particle ID][0]
-    FHist_(2,List<double**>(FHistSize_,NULL)),     // FHist_[k={1,2}-1][ndt in past][particle ID][dim]
-    gH0_(NULL),
-    tRef_(NULL),
-    mRef_(NULL),
-    lRef_(NULL),
+    nInt_(propsDict_.lookupOrDefault<int>("nIntegral", -1)),
+    discOrder_(propsDict_.lookupOrDefault<int>("discretisationOrder", 1)),
+    ddtUrelHistSize_(nInt_+discOrder_),
+    rHistSize_(nInt_),
+    FHistSize_(2*discOrder_),
+    ddtUrelHistRegName_(typeName + "ddtUrelHist"),  // indexed as: ddtUrelHist_[particle ID][iDim*ddtUrelHistSize_+iHist]
+    rHistRegName_(typeName + "rHist"),              // indexed as: rHist_[particle ID][iHist]
+    FHistRegName_(typeName + "FHist"),              // indexed as: ddtUrelHist_[particle ID][iDim*FHistSize_*2+iHist*2+k]
+    gH0RegName_(typeName + "gH0"),
+    tRefRegName_(typeName + "tRef"),
+    mRefRegName_(typeName + "mRef"),
+    lRefRegName_(typeName + "lRef"),
     Urel_
     (   IOobject
         (
@@ -115,12 +115,12 @@ ParmarBassetForce::ParmarBassetForce
         )
     )
 {
-
     // init force sub model
     setForceSubModels(propsDict_);
     // define switches which can be read from dict
     forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
     forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
     forceSubM(0).readSwitches();
 
     //Extra switches/settings
@@ -129,15 +129,23 @@ ParmarBassetForce::ParmarBassetForce
     if (discOrder_ < 1 || discOrder_ > 2)
         FatalError << "Parmar Basset Force: Discretisation order > 2 not implemented!" << abort(FatalError);
 
+    if (nInt_ < 1)
+        FatalError << "Parmar Basset Force: nIntegral missing or invalid, must be > 1" << abort(FatalError);
+
+    // allocate particle properties
+    particleCloud_.registerParticleProperty<double**>(ddtUrelHistRegName_,3*ddtUrelHistSize_,NOTONCPU,false);
+    particleCloud_.registerParticleProperty<double**>(      rHistRegName_,        rHistSize_,NOTONCPU,false);
+    particleCloud_.registerParticleProperty<double**>(      FHistRegName_,      6*FHistSize_,NOTONCPU,false);
+    particleCloud_.registerParticleProperty<double**>(        gH0RegName_,                1 ,NOTONCPU,false);
+    particleCloud_.registerParticleProperty<double**>(       tRefRegName_,                1 ,NOTONCPU,false);
+    particleCloud_.registerParticleProperty<double**>(       mRefRegName_,                1 ,NOTONCPU,false);
+    particleCloud_.registerParticleProperty<double**>(       lRefRegName_,                1 ,NOTONCPU,false);
+
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("ParmarBassetForce"); //first entry must the be the force
     particleCloud_.probeM().vectorFields_.append("Urel");
     particleCloud_.probeM().vectorFields_.append("ddtUrel");
-    //
-    particleCloud_.probeM().vectorFields_.append("UrelNoSmooth");
-    particleCloud_.probeM().vectorFields_.append("ddtUrelNoSmooth");
-    //
     particleCloud_.probeM().vectorFields_.append("Fshort");
     particleCloud_.probeM().vectorFields_.append("Flong1");
     particleCloud_.probeM().vectorFields_.append("Flong2");
@@ -158,20 +166,6 @@ ParmarBassetForce::ParmarBassetForce
 
 ParmarBassetForce::~ParmarBassetForce()
 {
-    particleCloud_.dataExchangeM().destroy(gH0_,    1);
-    particleCloud_.dataExchangeM().destroy(tRef_,   1);
-    particleCloud_.dataExchangeM().destroy(mRef_,   1);
-    particleCloud_.dataExchangeM().destroy(lRef_,   1);
-
-    for (int i=0; i<nHist_; i++)
-    {
-        particleCloud_.dataExchangeM().destroy(ddtUrelHist_[i],3);
-        particleCloud_.dataExchangeM().destroy(rHist_      [i],1);
-    }
-
-    for (int i=0; i<FHistSize_; i++)
-        for (int k=0; k<2; k++)
-            particleCloud_.dataExchangeM().destroy(FHist_[k][i],3);
 
 }
 
@@ -180,16 +174,18 @@ ParmarBassetForce::~ParmarBassetForce()
 
 void ParmarBassetForce::setForce() const
 {
-    // allocate arrays
-    if(particleCloud_.numberOfParticlesChanged())
-        reAllocArrays();
+    double**& ddtUrelHist_ = particleCloud_.getParticlePropertyRef<double**>(ddtUrelHistRegName_);
+    double**& rHist_       = particleCloud_.getParticlePropertyRef<double**>(      rHistRegName_);
+    double**& FHist_       = particleCloud_.getParticlePropertyRef<double**>(      FHistRegName_);
+    double**& gH0_         = particleCloud_.getParticlePropertyRef<double**>(        gH0RegName_);
+    double**& tRef_        = particleCloud_.getParticlePropertyRef<double**>(       tRefRegName_);
+    double**& mRef_        = particleCloud_.getParticlePropertyRef<double**>(       mRefRegName_);
+    double**& lRef_        = particleCloud_.getParticlePropertyRef<double**>(       lRefRegName_);
 
     vector position(0,0,0);
     vector Urel(0,0,0);
     vector ddtUrel(0,0,0);
 
-    scalar t0min    = 0.00;
-
     const volScalarField& nufField = forceSubM(0).nuField();
     const volScalarField& rhoField = forceSubM(0).rhoField();
 
@@ -220,26 +216,17 @@ void ParmarBassetForce::setForce() const
     Urel_ = Us_ - U_;
     ddtUrel_ = fvc::ddt(Us_) - fvc::ddt(U_) - (Us_ & fvc::grad(U_));
 
-    //
-    volVectorField UrelNoSmooth_ = Urel_;
-    volVectorField ddtUrelNoSmooth_ = ddtUrel_;
-    //
-
     smoothingM().smoothen(Urel_);
     smoothingM().smoothen(ddtUrel_);
 
     interpolationCellPoint<vector> UrelInterpolator_(Urel_);
     interpolationCellPoint<vector> ddtUrelInterpolator_(ddtUrel_);
 
-    //
-    interpolationCellPoint<vector> UrelNoSmoothInterpolator_(UrelNoSmooth_);
-    interpolationCellPoint<vector> ddtUrelNoSmoothInterpolator_(ddtUrelNoSmooth_);
-    //
-
     for(int index = 0;index <  particleCloud_.numberOfParticles(); index++)
     {
             vector ParmarBassetForce(0,0,0);
             vector Fshort(0,0,0);
+            vector Flong[2]={vector::zero, vector::zero};
             label  cellI = particleCloud_.cellIDs()[index][0];
 
             if (cellI > -1) // particle Found
@@ -268,12 +255,13 @@ void ParmarBassetForce::setForce() const
 
                 if (gH0_[index][0]!=NOTONCPU)
                 {
-                    r = pow(gH0_[index][0]/gH,1.5); // Eq. 3.4
+                    scalar gHratio = gH0_[index][0]/gH;
+                    r = gHratio*sqrt(gHratio); // gHratio^1.5, Eq. 3.4
 
                     if (r<0.25 || r>2.0)
                     {
-                        gH0_[index][0]            = NOTONCPU; //reset reference
-                        ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
+                        gH0_[index][0]         = NOTONCPU; //reset reference
+                        ddtUrelHist_[index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
                     }
 
                 }
@@ -285,7 +273,7 @@ void ParmarBassetForce::setForce() const
                     scalar Vs   = rs*rs*rs*M_PI*4/3;
                     scalar mRef = Vs*rhoField[cellI] * gH * 5.2863;             // 9/(2*sqrt(pi))*(256/pi)^(1/6) = 5.2863 (Eq. 3.2)
 
-                    gH0_[index][0] = gH;
+                    gH0_[index][0]  = gH;
                     tRef_[index][0] = tRef;
                     mRef_[index][0] = mRef;
                     lRef_[index][0] = rs;
@@ -298,29 +286,15 @@ void ParmarBassetForce::setForce() const
                 scalar dt     = U_.mesh().time().deltaT().value() / tRef_[index][0];  // dim.less
                 scalar t0     = nInt_*dt; // dim.less
 
-                //********* update histories *********//
-
                 // non-dimensionlise
                 Urel    /= mps;
                 ddtUrel /= mpss;
 
-                // update ddtUrel  and r history
-                update_ddtUrelHist(ddtUrel,index); // add current dim.less ddtUrel to history
-                update_rHist(r,index); // add current r to history
-
-                // warning and reset for too small t0
-                if (t0<t0min)
-                {
-                    Pout << "ParmarBassetForce WARNING: t0 = " << t0 << " at ID = " << index <<  endl;
-                    gH0_[index][0]            = NOTONCPU; //reset reference
-                    ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
-                }
-
                 // check length of known history
-                int nKnown = 0;
-                for (int j=0; j<nHist_; j++) // loop over past times
+                int nKnown = 1; // we always know the current step
+                for (int j=0; j<ddtUrelHistSize_; j++) // loop over past times
                 {
-                    if (ddtUrelHist_[j][index][0] == NOTONCPU)
+                    if (ddtUrelHist_[index][j] == NOTONCPU)
                         break;
                     else
                         nKnown++;
@@ -332,117 +306,92 @@ void ParmarBassetForce::setForce() const
                 int nShort = min(nKnown,nInt_+1);
 
                 // int_0^dt K(r,xi) dxi * ddtU(t) dxi (singularity treated by assuming constant acceleration)
-                if (nShort>0)
-                {
-                    for (int i=0; i<3; i++) // loop over dimensions
-                        Fshort[i] = -calculateK0(r,dt) * ddtUrelHist_[0][index][i];
-                }
+                for (int i=0; i<3; i++) // loop over dimensions
+                    Fshort[i] = -calculateK0(r,dt) * ddtUrel[i];
 
                 // int_dt^t0 K(r,xi) * ddtU(t-xi) dxi (trapezoid rule)
                 if (nShort>2)
                 {
-                    for (int j=1; j<nShort; j++)
+                    for (int j=0; j<(nShort-1); j++) // we don't use the current step here, hence nShort-1
                     {
-                        scalar xi = j*dt;
-                        scalar K  = pow((pow(xi,.25) + rHist_[j][index][0]*xi),-2.); // Eq. 3.4
+                        scalar xi = (j+1)*dt;
+                        scalar invsqrtK = sqrt(sqrt(xi)) + rHist_[index][j]*xi; // K^-0.5
+                        scalar K  = 1./(invsqrtK*invsqrtK); // Eq. 3.4
 
                         for (int i=0; i<3; i++) // loop over dimensions
-                            Fshort[i] -= trapWeight(j,nShort) * K * ddtUrelHist_[j][index][i] * dt;
+                            Fshort[i] -= trapWeight(j,nShort-1) * K * ddtUrelHist_[index][i*ddtUrelHistSize_+j] * dt;
                     }
                 }
 
                 //********* long term force computing (differential form) *********//
 
-                // update F1, F2 history
-                update_FHist(vector::zero,vector::zero,index);
-
                 // initialise ddtUrel(t0) and Flong(:) as 0 and r(t0) as 1
                 if (nKnown == nInt_)
                 {
-                    for (int j=nInt_; j<nHist_; j++) // loop over past times
-                    {
-                        rHist_[j][index][0] = 1.;
+                    // initialise the histories beyond nInt
+                    for (int j=nInt_-1; j<rHistSize_; j++) // loop over past times
+                        rHist_[index][j] = 1.;
+
+                    for (int j=nInt_-1; j<ddtUrelHistSize_; j++) // loop over past times
                         for (int i=0; i<3; i++) // loop over dimensions
-                            ddtUrelHist_[j][index][i] = 0.0;
-                    }
+                            ddtUrelHist_[index][i*ddtUrelHistSize_+j] = 0.0;
 
                     for (int k=0; k<2; k++) // loop over F1, F2
                         for (int j=0; j<FHistSize_; j++) // loop over past times
                             for (int i=0; i<3; i++) // loop over dimensions
-                                FHist_[k][j][index][i] = 0.0;
-                    nKnown = nHist_;
+                                FHist_[index][i*FHistSize_*2+j*2+k] = 0.0;
+                    nKnown = ddtUrelHistSize_+1;
                 }
 
                 // solve ODEs
-                if (nKnown == nHist_)
+                if (nKnown == ddtUrelHistSize_+1)
                 {
                     for (int k=0; k<2; k++) // loop over F1, F2
                     {
                         //calculate coefficients
                         double C[4];
-                        calculateCoeffs(k,t0,rHist_[nInt_][index][0],c,chi,C);
+                        calculateCoeffs(k,t0,rHist_[index][nInt_-1],c,chi,C);
 
                         // solve Eq. 3.20
-                        solveFlongODE(k,C,dt,index);
+                        Flong[k] = solveFlongODE(FHist_,ddtUrelHist_,k,C,dt,index);
                     }
                 }
 
+                //********* update histories *********//
+                update_ddtUrelHist(ddtUrelHist_,ddtUrel,index); // add current dim.less ddtUrel to history
+                update_rHist(rHist_,r,index); // add current r to history
+                update_FHist(FHist_,Flong[0],Flong[1],index);
+
                 //********* total force *********//
 
                 // sum and convert to N
-                for (int i=0; i<3; i++) // loop over dimensions
-                {
-                    ParmarBassetForce[i] = Fshort[i];
-                    for (int k=0; k<2; k++) // loop over F1, F2
-                        ParmarBassetForce[i] += FHist_[k][0][index][i];
-                }
+                ParmarBassetForce = Fshort;
+                for (int k=0; k<2; k++) // loop over F1, F2
+                    ParmarBassetForce += Flong[k];
                 ParmarBassetForce *= newton;
 
+                if (forceSubM(0).verbose() && index >= 0 && index < 2)
+                {
+                    Pout << "cellI = " << cellI << endl;
+                    Pout << "index = " << index << endl;
+                    Pout << "Fshort = " << Fshort*newton << endl;
+                    Pout << "Flong1 = " << Flong[0]*newton << endl;
+                    Pout << "Flong2 = " << Flong[1]*newton << endl;
+                    Pout << "Ftotal = " << ParmarBassetForce << endl;
+                }
+
                 // Set value fields and write the probe
                 if(probeIt_)
                 {
                     scalar ReRef = 0.75/(gH0_[index][0]-0.105);
 
-                    vector Flong1(0,0,0);
-                    vector Flong2(0,0,0);
-
-                    for (int i=0; i<3; i++) // loop over dimensions
-                    {
-                        Flong1[i] = FHist_[0][0][index][i];
-                        Flong2[i] = FHist_[1][0][index][i];
-                    }
-
-                    //
-                    // relative velocity (m/s)
-                    vector UrelNoSmooth;
-                    vector ddtUrelNoSmooth;
-
-                    if(forceSubM(0).interpolation())
-                        UrelNoSmooth = UrelNoSmoothInterpolator_.interpolate(position,cellI);
-                    else
-                        UrelNoSmooth = UrelNoSmooth_[cellI];
-
-                    // acceleration (m/s2)
-                    if(forceSubM(0).interpolation())
-                        ddtUrelNoSmooth = ddtUrelNoSmoothInterpolator_.interpolate(position,cellI);
-                    else
-                        ddtUrelNoSmooth = ddtUrelNoSmooth_[cellI];
-
-                    UrelNoSmooth /= mps;
-                    ddtUrelNoSmooth /= mpss;
-                    //
-
                     #include "setupProbeModelfields.H"
                     vValues.append(ParmarBassetForce);           //first entry must the be the force
                     vValues.append(Urel);
                     vValues.append(ddtUrel);
-                    //
-                    vValues.append(UrelNoSmooth);
-                    vValues.append(ddtUrelNoSmooth);
-                    //
                     vValues.append(Fshort);
-                    vValues.append(Flong1);
-                    vValues.append(Flong2);
+                    vValues.append(Flong[0]);
+                    vValues.append(Flong[1]);
                     sValues.append(ReRef);
                     sValues.append(tRef_[index][0]);
                     sValues.append(mRef_[index][0]);
@@ -456,8 +405,8 @@ void ParmarBassetForce::setForce() const
             else
             {
                 // not on CPU
-                gH0_[index][0]            = NOTONCPU; //reset reference
-                ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
+                gH0_[index][0]         = NOTONCPU; //reset reference
+                ddtUrelHist_[index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
             }
 
             // write particle based data to global array
@@ -465,32 +414,11 @@ void ParmarBassetForce::setForce() const
     }
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::ParmarBassetForce::reAllocArrays() const
-{
-    Pout << "ParmarBassetForce::reAllocArrays..." << endl;
-
-    particleCloud_.dataExchangeM().allocateArray(gH0_,     NOTONCPU,1);
-    particleCloud_.dataExchangeM().allocateArray(tRef_,    NOTONCPU,1);
-    particleCloud_.dataExchangeM().allocateArray(mRef_,    NOTONCPU,1);
-    particleCloud_.dataExchangeM().allocateArray(lRef_,    NOTONCPU,1);
-
-    for (int i=0; i<nHist_; i++)
-    {
-        particleCloud_.dataExchangeM().allocateArray(ddtUrelHist_[i],NOTONCPU,3);
-        particleCloud_.dataExchangeM().allocateArray(rHist_      [i],NOTONCPU,1);
-    }
-
-    for (int i=0; i<2*discOrder_+1; i++)
-        for (int k=0; k<2; k++)
-            particleCloud_.dataExchangeM().allocateArray(FHist_[k][i],0.0,3);
-}
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 scalar Foam::ParmarBassetForce::calculateK0(scalar r, scalar dt) const
 {
     scalar cbrtr = cbrt(r); // cube root of r
-    scalar gamma = cbrtr*pow(dt,0.25);
+    scalar gamma = cbrtr*sqrt(sqrt(dt));
 
     /*
     scalar K0 = 2./(9.*pow(r,0.666)) *
@@ -513,44 +441,44 @@ scalar Foam::ParmarBassetForce::calculateK0(scalar r, scalar dt) const
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 scalar Foam::ParmarBassetForce::trapWeight(int i, int n) const
 {
-    if ( (i==1) || (i==(n-1)) )
+    if ( (i==0) || (i==(n-1)) )
         return 0.5;
     else
         return 1.0;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::ParmarBassetForce::update_ddtUrelHist(const vector& ddtUrel, int index) const
+void Foam::ParmarBassetForce::update_ddtUrelHist(double**& ddtUrelHist_, const vector& ddtUrel, int index) const
 {
     for (int i=0; i<3; i++) // loop over dimensions
     {
-        for (int j=nHist_-1; j>0; j--) // loop over past times
-            ddtUrelHist_[j][index][i] = ddtUrelHist_[j-1][index][i];
+        for (int j=ddtUrelHistSize_-1; j>0; j--) // loop over past times
+            ddtUrelHist_[index][i*ddtUrelHistSize_+j] = ddtUrelHist_[index][i*ddtUrelHistSize_+j-1];
 
-        ddtUrelHist_[0][index][i] = ddtUrel[i];
+        ddtUrelHist_[index][i*ddtUrelHistSize_] = ddtUrel[i];
     }
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::ParmarBassetForce::update_rHist(scalar r, int index) const
+void Foam::ParmarBassetForce::update_rHist(double**& rHist_, scalar r, int index) const
 {
-    for (int j=nHist_-1; j>0; j--) // loop over past times
-        rHist_[j][index][0] = rHist_[j-1][index][0];
+    for (int j=rHistSize_-1; j>0; j--) // loop over past times
+        rHist_[index][j] = rHist_[index][j-1];
 
-    rHist_[0][index][0] = r;
+    rHist_[index][0] = r;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::ParmarBassetForce::update_FHist(const vector& F1, const vector& F2, int index) const
+void Foam::ParmarBassetForce::update_FHist(double**& FHist_, const vector& F1, const vector& F2, int index) const
 {
     for (int i=0; i<3; i++) // loop over dimensions
     {
         for (int k=0; k<2; k++) // loop over F1, F2
             for (int j=FHistSize_-1; j>0; j--) // loop over past times
-                FHist_[k][j][index][i] = FHist_[k][j-1][index][i];
+                FHist_[index][i*FHistSize_*2+j*2+k] = FHist_[index][i*FHistSize_*2+(j-1)*2+k];
 
-        FHist_[0][0][index][i] = F1[i];
-        FHist_[1][0][index][i] = F2[i];
+        FHist_[index][i*FHistSize_*2  ] = F1[i];
+        FHist_[index][i*FHistSize_*2+1] = F2[i];
     }
 }
 
@@ -564,18 +492,20 @@ void Foam::ParmarBassetForce::calculateCoeffs(int k, scalar t0, scalar r, double
     }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::ParmarBassetForce::solveFlongODE(int k, double C[4], scalar dt, int index) const
+vector Foam::ParmarBassetForce::solveFlongODE(double**& FHist_, double**& ddtUrelHist_, int k, double C[4], scalar dt, int index) const
 {
+    vector Flong = vector::zero;
+
     if (discOrder_==1)
     {
         for (int i=0; i<3; i++) // loop over dimensions
         {
-            FHist_[k][0][index][i] =
+            Flong[i] = 
                     (
-                        (     C[1]/dt+2/(dt*dt)) * FHist_[k][1][index][i]
-                      - (             1/(dt*dt)) * FHist_[k][2][index][i]
-                      - (C[2]+C[3]/dt        ) * ddtUrelHist_[nInt_  ][index][i]
-                      + (     C[3]/dt        ) * ddtUrelHist_[nInt_+1][index][i]
+                        (     C[1]/dt+2/(dt*dt)) * FHist_[index][i*FHistSize_*2  +k]
+                      - (             1/(dt*dt)) * FHist_[index][i*FHistSize_*2+2+k]
+                      - (C[2]+C[3]/dt        ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_-1]
+                      + (     C[3]/dt        ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_  ]
                     ) / (C[0]+C[1]/dt+1/(dt*dt)); // Eq. 3.20 using first order temporal discretisation
         }
     }
@@ -583,39 +513,20 @@ void Foam::ParmarBassetForce::solveFlongODE(int k, double C[4], scalar dt, int i
     {
         for (int i=0; i<3; i++) // loop over dimensions
         {
-            FHist_[k][0][index][i] =
+            Flong[i] = 
                     (
-                        (       4*C[1]/(2*dt) + 24/(4*dt*dt)) * FHist_[k][1][index][i]
-                      - (         C[1]/(2*dt) + 22/(4*dt*dt)) * FHist_[k][2][index][i]
-                      + (                        8/(4*dt*dt)) * FHist_[k][3][index][i]
-                      - (                        1/(4*dt*dt)) * FHist_[k][4][index][i]
+                        (       4*C[1]/(2*dt) + 24/(4*dt*dt)) * FHist_[index][i*FHistSize_*2  +k]
+                      - (         C[1]/(2*dt) + 22/(4*dt*dt)) * FHist_[index][i*FHistSize_*2+2+k]
+                      + (                        8/(4*dt*dt)) * FHist_[index][i*FHistSize_*2+4+k]
+                      - (                        1/(4*dt*dt)) * FHist_[index][i*FHistSize_*2+6+k]
 
-                      - (C[2] + 3*C[3]/(2*dt)               ) * ddtUrelHist_[nInt_  ][index][i]
-                      + (       4*C[3]/(2*dt)               ) * ddtUrelHist_[nInt_+1][index][i]
-                      - (         C[3]/(2*dt)               ) * ddtUrelHist_[nInt_+2][index][i]
+                      - (C[2] + 3*C[3]/(2*dt)               ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_-1]
+                      + (       4*C[3]/(2*dt)               ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_  ]
+                      - (         C[3]/(2*dt)               ) * ddtUrelHist_[index][i*ddtUrelHistSize_+nInt_+1]
                     ) / (C[0] + 3*C[1]/(2*dt) +  9/(4*dt*dt)); // Eq. 3.20 using second order temporal discretisation
         }
     }
-}
-
-void Foam::ParmarBassetForce::rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const
-{
-    for (int i=0; i<3; i++) // loop over dimensions
-    {
-        // rescale ddtU history
-        for (int j=0; j<nHist_; j++) // loop over past times
-            if (ddtUrelHist_[j][index][i] != NOTONCPU)
-                ddtUrelHist_[j][index][i] *= lScale/(tScale*tScale);
-
-        // rescale F1, F2 history
-        for (int k=0; k<2; k++) // loop over F1, F2
-            for (int j=0; j<FHistSize_; j++) // loop over past times
-                FHist_[k][j][index][i] *= mScale*lScale/(tScale*tScale);
-    }
-    // rescale r history
-    for (int j=0; j<nHist_; j++) // loop over past times
-        if (rHist_[j][index][0] != NOTONCPU)
-            rHist_[j][index][0] /= pow(rScale,1.5);
+    return Flong;
 }
 
 } // End namespace Foam

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

@@ -59,35 +59,37 @@ class ParmarBassetForce
 private:
     dictionary propsDict_;
 
-    word velFieldName_;
+    const word velFieldName_;
 
     const volVectorField& U_;
 
-    word UsFieldName_;
+    const word UsFieldName_;
 
     const volVectorField& Us_;
 
-    label nInt_;                            //no of timesteps to solve full integral
+    const int nInt_;                            //no of timesteps to solve full integral
 
-    label discOrder_;                       //ODE discretisation order
+    const int discOrder_;                       //ODE discretisation order
 
-    label nHist_;                           //no of timesteps to save ddtUrel history for
+    const int ddtUrelHistSize_;                 //no of timesteps to save ddtUrel history for
 
-    label FHistSize_;
+    const int rHistSize_;                       //no of timesteps to save r history for
 
-    mutable List<double**> ddtUrelHist_;
+    const int FHistSize_;                       //no of timesteps to save Flong history for
 
-    mutable List<double**> rHist_;
+    const word ddtUrelHistRegName_;
 
-    mutable List<List<double**>> FHist_;
+    const word rHistRegName_;
 
-    mutable double** gH0_;
+    const word FHistRegName_;
 
-    mutable double** tRef_;
+    const word gH0RegName_;
 
-    mutable double** mRef_;
+    const word tRefRegName_;
 
-    mutable double** lRef_;
+    const word mRefRegName_;
+
+    const word lRefRegName_;
 
     mutable volVectorField Urel_;
 
@@ -101,17 +103,15 @@ private:
 
     scalar trapWeight(int i, int n) const;
 
-    void update_ddtUrelHist(const vector& ddtUrel, int index) const;
+    void update_ddtUrelHist(double**& ddtUrelHist_, const vector& ddtUrel, int index) const;
 
-    void update_rHist(scalar r, int index) const;
+    void update_rHist(double**& rHist_, scalar r, int index) const;
 
-    void update_FHist(const vector& F1, const vector& F2, int index) const;
+    void update_FHist(double**& FHist_, const vector& F1, const vector& F2, int index) const;
 
     void calculateCoeffs(int k, scalar t0, scalar r, double c[2][4][3], double chi[2][4][2], double C[4]) const;
 
-    void solveFlongODE(int k, double C[4], scalar dt, int index) const;
-
-    void rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const;
+    vector solveFlongODE(double**& FHist_, double**& ddtUrelHist_, int k, double C[4], scalar dt, int index) const;
 
 public:
 
@@ -136,9 +136,6 @@ public:
     // Member Functions
         void setForce() const;
 
-
-        void reAllocArrays() const;
-
         inline const smoothingModel& smoothingM() const
         {
             return smoothingModel_;

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

@@ -5,7 +5,8 @@
     www.cfdem.com
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
-                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright 2012-2015 DCS Computing GmbH, Linz
+                                Copyright 2015-     JKU Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -66,10 +67,14 @@ ShirgaonkarIB::ShirgaonkarIB
     verbose_(propsDict_.found("verbose")),
     twoDimensional_(propsDict_.found("twoDimensional")),
     depth_(1),
+    multisphere_(propsDict_.found("multisphere")), //  drag for a multisphere particle
+    useTorque_(propsDict_.found("useTorque")),
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     pressureFieldName_(propsDict_.lookup("pressureFieldName")),
-    p_(sm.mesh().lookupObject<volScalarField> (pressureFieldName_))
+    p_(sm.mesh().lookupObject<volScalarField> (pressureFieldName_)),
+    solidVolFractionName_(multisphere_?propsDict_.lookup("solidVolFractionName"):word::null),
+    lambda_(multisphere_?sm.mesh().lookupObject<volScalarField>(solidVolFractionName_):volScalarField::null())
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
@@ -109,6 +114,7 @@ void ShirgaonkarIB::setForce() const
 
     label cellI;
     vector drag;
+    vector torque;
 
     volVectorField h=forceSubM(0).IBDragPerV(U_,p_);
 
@@ -119,6 +125,7 @@ void ShirgaonkarIB::setForce() const
         //if(mask[index][0])
         //{
             drag=vector::zero;
+            torque=vector::zero;
 
             for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
             {
@@ -128,6 +135,12 @@ void ShirgaonkarIB::setForce() const
                 if (cellI > -1) // particle Found
                 {
                     drag += h[cellI]*h.mesh().V()[cellI];
+                    if(useTorque_)
+                    {
+                        vector rc = particleCloud_.mesh().C()[cellI];
+                        vector positionCenter = particleCloud_.position(index);
+                        torque += (rc - positionCenter)^h[cellI]*h.mesh().V()[cellI];
+                    }
                 }
 
             }
@@ -150,10 +163,38 @@ void ShirgaonkarIB::setForce() const
                               << ","            << particleCloud_.impForces()[index][1]
                               << ","            << particleCloud_.impForces()[index][2]
                               << endl;
+
+            if(useTorque_)
+            {
+                for(int j=0;j<3;j++) particleCloud_.DEMTorques()[index][j] = torque[j]; // Adding to the particle torque;
+                if(verbose_) Info << "DEMTorques = " << particleCloud_.DEMTorques()[index][0] << "," << particleCloud_.DEMTorques()[index][1] << "," << particleCloud_.DEMTorques()[index][2] << endl;
+            }
         //}
     }
+
+    // Calculate the force if the particle is multisphere template
+    if(multisphere_) calcForce();
 }
 
+void ShirgaonkarIB::calcForce() const
+{
+    vector dragMS;
+
+    volVectorField h=forceSubM(0).IBDragPerV(U_,p_);
+
+    dragMS = vector::zero;
+
+    forAll(lambda_,cellI)
+    {
+        if(lambda_[cellI] > 0) dragMS += h[cellI]*h.mesh().V()[cellI];
+        else dragMS = dragMS;
+    }
+
+    Pout << "Drag force on particle clump = " << dragMS[0] << ", " << dragMS[1] << ", " << dragMS[2] << endl;
+}
+
+
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -5,7 +5,8 @@
     www.cfdem.com
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
-                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright 2012-2015 DCS Computing GmbH, Linz
+                                Copyright 2015-     JKU Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -66,6 +67,10 @@ private:
 
     bool depth_;
 
+    const bool multisphere_;
+
+    const bool useTorque_;
+
     word velFieldName_;
 
     const volVectorField& U_;
@@ -74,6 +79,12 @@ private:
 
     const volScalarField& p_;
 
+    word solidVolFractionName_;
+
+    const volScalarField& lambda_;
+
+    void calcForce() const;
+
 public:
 
     //- Runtime type information

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

@@ -61,6 +61,7 @@ gradPForceSmooth::gradPForceSmooth
     U_(sm.mesh().lookupObject<volVectorField> (velocityFieldName_)),
     useRho_(false),
     useU_(false),
+    temporalSmoothing_(false),
     addedMassCoeff_(0.0),
     smoothingModel_
     (
@@ -96,7 +97,7 @@ gradPForceSmooth::gradPForceSmooth
 
     if (modelType_ == "B")
     {
-        FatalError <<"using  model gradPForceSmooth with model type B is not valid\n" << abort(FatalError);
+        FatalError <<"using model gradPForceSmooth with model type B is not valid\n" << abort(FatalError);
     }
     else if (modelType_ == "Bfull")
     {
@@ -135,6 +136,11 @@ gradPForceSmooth::gradPForceSmooth
     particleCloud_.probeM().scalarFields_.append("Vs");
     particleCloud_.probeM().scalarFields_.append("rho");
     particleCloud_.probeM().writeHeader();
+
+    if ((smoothingM().type() == "temporalSmoothing") || (smoothingM().type() == "constDiffAndTemporalSmoothing"))
+    {
+        temporalSmoothing_ = true;
+    }
 }
 
 
@@ -150,6 +156,11 @@ gradPForceSmooth::~gradPForceSmooth()
 void gradPForceSmooth::setForce() const
 {
     volVectorField gradPField = fvc::grad(p_);
+
+    // if temporal smoothing is used, let pSmooth evolve on its own - smoothing model itself looks up p
+    // else, set pSmooth to current p and apply spatial smoothing
+    if (!temporalSmoothing_) pSmooth_ = p_;
+
     if (pFieldName_ == "p_rgh")
     {
         const volScalarField& rho_ = particleCloud_.mesh().lookupObject<volScalarField>("rho");

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

@@ -62,6 +62,8 @@ private:
 
     bool useU_;          // if false: substitution p=0.5*rho*U^2
 
+    bool temporalSmoothing_;
+
     mutable double addedMassCoeff_; //added mass coefficient
 
     autoPtr<smoothingModel> smoothingModel_;

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

@@ -103,7 +103,6 @@ interface::~interface()
 
 void interface::setForce() const
 {
-Info << "interface::setForce" << endl;
     for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
     {
         //if(mask[index][0])
@@ -193,7 +192,6 @@ Info << "interface::setForce" << endl;
             } // end if particle found on proc domain
         //}// end if in mask
     }// end loop particles
-Info << "interface::setForce - done" << endl;
 }
 
 

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

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+Description
+    transfer fluid properties to LIGGGHTS
+
+SourceFiles
+    transferFluidProperties.C
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "transferFluidProperties.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(transferFluidProperties, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    transferFluidProperties,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+transferFluidProperties::transferFluidProperties
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props"))
+{
+    particleCloud_.registerParticleProperty<double**>("fluidDensity",1);
+    particleCloud_.registerParticleProperty<double**>("fluidViscosity",1);
+
+    // init force sub model
+    setForceSubModels(propsDict_);
+    // define switches which can be read from dict
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).readSwitches();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+transferFluidProperties::~transferFluidProperties()
+{
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
+
+void transferFluidProperties::setForce() const
+{
+    double**& fluidDensity_   = particleCloud_.getParticlePropertyRef<double**>("fluidDensity");
+    double**& fluidViscosity_ = particleCloud_.getParticlePropertyRef<double**>("fluidViscosity");
+
+    const volScalarField& rhoField = forceSubM(0).rhoField();
+    const volScalarField& nufField = forceSubM(0).nuField();
+
+    interpolationCellPoint<scalar> rhoInterpolator_(rhoField);
+    interpolationCellPoint<scalar> nufInterpolator_(nufField);
+
+    label cellI = 0;
+    double rho = 0.;
+    double nuf = 0.;
+    vector position(0,0,0);
+
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
+    {
+        cellI = particleCloud_.cellIDs()[index][0];
+        if (cellI >= 0)
+        {
+            if(forceSubM(0).interpolation())
+            {
+                position = particleCloud_.position(index);
+                rho = rhoInterpolator_.interpolate(position,cellI);
+                nuf = nufInterpolator_.interpolate(position,cellI);
+            }
+            else
+            {
+                rho = rhoField[cellI];
+                nuf = nufField[cellI];
+            }
+
+            fluidDensity_[index][0]   = rho;
+            fluidViscosity_[index][0] = nuf*rho;
+        }
+    }
+
+    particleCloud_.dataExchangeM().giveData("fluidDensity","scalar-atom",fluidDensity_);
+    particleCloud_.dataExchangeM().giveData("fluidViscosity","scalar-atom",fluidViscosity_);
+    
+    if (forceSubM(0).verbose()) Info << "give data done" << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

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

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

@@ -99,7 +99,12 @@ virtualMassForce::virtualMassForce
         )
     )
 {
-    particleCloud_.registerParticleProperty<double**>(UrelOldRegName_,3,NOTONCPU,false);
+    // allocate UrelOld only if needed
+    int UrelOldSize = 0;
+    if(!splitUrelCalculation_ || !useUs_ )
+        UrelOldSize = 3;
+
+    particleCloud_.registerParticleProperty<double**>(UrelOldRegName_,UrelOldSize,NOTONCPU,false);
 
     // init force sub model
     setForceSubModels(propsDict_);

src/lagrangian/cfdemParticle/subModels/massTransferModel/massTransferModel/massTransferModel.H

@@ -105,8 +105,6 @@ public:
 
         virtual void postFlow() {}
 
-        virtual void solve() {}
-
         const volScalarField& D0Field() const;
 
         const volScalarField& CsField() const;

src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.C

@@ -219,7 +219,7 @@ void particleProbe::writeProbe(int index, Field<scalar> sValues, Field<vector> v
         *sPtr << setprecision(writePrecision_);
         forAll(vValues, iter)
         {
-            // if(!probeDebug_ && iter>0) break;
+            if(!probeDebug_ && iter>0) break;
             *sPtr << vValues[iter][0] << "   ";
             *sPtr << vValues[iter][1] << "   ";
             *sPtr << vValues[iter][2] << "   ";

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

@@ -68,6 +68,7 @@ constDiffAndTemporalSmoothing::constDiffAndTemporalSmoothing
     smoothingLength_(dimensionedScalar("smoothingLength",dimensionSet(0,1,0,0,0,0,0), readScalar(propsDict_.lookup("smoothingLength")))),
     smoothingLengthReferenceField_(dimensionedScalar("smoothingLengthReferenceField",dimensionSet(0,1,0,0,0,0,0), readScalar(propsDict_.lookup("smoothingLength")))),
     DT_("DT", dimensionSet(0,2,-1,0,0), 0.),
+    refFieldName_(propsDict_.lookupOrDefault<word>("refField","")),
     gamma_(readScalar(propsDict_.lookup("smoothingStrength"))),
     correctBoundary_(propsDict_.lookupOrDefault<bool>("correctBoundary",false)),
     verbose_(false)
@@ -111,12 +112,43 @@ void constDiffAndTemporalSmoothing::smoothen(volScalarField& fieldSrc) const
     double deltaT = sSmoothField.mesh().time().deltaTValue();
     DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT;
 
-    // do smoothing
-    solve
-    (
-        fvm::ddt(sSmoothField)
-       -fvm::laplacian(DT_, sSmoothField)
-    );
+    if(refFieldName_.empty())
+    {
+        // do spatial smoothing
+        solve
+        (
+            fvm::ddt(sSmoothField)
+           -fvm::laplacian(DT_, sSmoothField)
+        );
+
+        // create fields for temporal smoothing
+        volScalarField refField = sSmoothField;
+
+        sSmoothField.ref()=fieldSrc.oldTime().internalField();
+        sSmoothField.correctBoundaryConditions();
+        sSmoothField.oldTime()=fieldSrc.oldTime();
+        sSmoothField.oldTime().correctBoundaryConditions();
+
+        // do temporal smoothing
+        solve
+        (
+            fvm::ddt(sSmoothField)
+            -
+            gamma_/deltaT * (refField - fvm::Sp(1.0,sSmoothField))
+        );
+    }
+    else
+    {
+        const volScalarField& refField = particleCloud_.mesh().lookupObject<volScalarField>(refFieldName_);
+
+        // do spatial and temporal smoothing
+        solve
+        (
+            fvm::ddt(sSmoothField)
+           -fvm::laplacian(DT_, sSmoothField)
+           -gamma_/deltaT * (refField - fvm::Sp(1.0,sSmoothField))
+        );
+    }
 
     // bound sSmoothField_
     forAll(sSmoothField,cellI)
@@ -152,28 +184,43 @@ void constDiffAndTemporalSmoothing::smoothen(volVectorField& fieldSrc) const
     dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
     DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT.value();
 
-    // do spacial smoothing
-    solve
-    (
-        fvm::ddt(vSmoothField)
-       -fvm::laplacian(DT_, vSmoothField)
-    );
+    if(refFieldName_.empty())
+    {
+        // do spatial smoothing
+        solve
+        (
+            fvm::ddt(vSmoothField)
+           -fvm::laplacian(DT_, vSmoothField)
+        );
 
-    // create fields for temporal smoothing
-    volVectorField refField = vSmoothField;
+        // create fields for temporal smoothing
+        volVectorField refField = vSmoothField;
 
-    vSmoothField.ref()=fieldSrc.oldTime().internalField();
-    vSmoothField.correctBoundaryConditions();
-    vSmoothField.oldTime()=fieldSrc.oldTime();
-    vSmoothField.oldTime().correctBoundaryConditions();
+        vSmoothField.ref()=fieldSrc.oldTime().internalField();
+        vSmoothField.correctBoundaryConditions();
+        vSmoothField.oldTime()=fieldSrc.oldTime();
+        vSmoothField.oldTime().correctBoundaryConditions();
 
-    // do temporal smoothing
-    solve
-    (
-        fvm::ddt(vSmoothField)
-        -
-        gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
-    );
+        // do temporal smoothing
+        solve
+        (
+            fvm::ddt(vSmoothField)
+            -
+            gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
+        );
+    }
+    else
+    {
+        const volVectorField& refField = particleCloud_.mesh().lookupObject<volVectorField>(refFieldName_);
+
+        // do spatial and temporal smoothing
+        solve
+        (
+            fvm::ddt(vSmoothField)
+           -fvm::laplacian(DT_, vSmoothField)
+           -gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
+        );
+    }
 
     // create temporally smoothened boundary field
     if (correctBoundary_)

src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.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
@@ -65,6 +65,7 @@ private:
     dimensionedScalar smoothingLength_;
     dimensionedScalar smoothingLengthReferenceField_;
     mutable dimensionedScalar DT_;
+    word refFieldName_;
     scalar gamma_;
     bool correctBoundary_;
     bool verbose_;
@@ -93,7 +94,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

@@ -130,7 +130,7 @@ void Foam::temporalSmoothing::smoothen(volVectorField& fieldSrc) const
     vSmoothField.oldTime()=fieldSrc;
     vSmoothField.oldTime().correctBoundaryConditions();
 
-    volVectorField refField = particleCloud_.mesh().lookupObject<volVectorField>(refFieldName_);
+    const volVectorField& refField = particleCloud_.mesh().lookupObject<volVectorField>(refFieldName_);
 
     dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
     solve

src/lagrangian/cfdemParticleComp/Make/files

@@ -90,6 +90,7 @@ $(forceModels)/directedDiffusiveRelaxation/directedDiffusiveRelaxation.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
 $(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
+$(forceModels)/transferFluidProperties/transferFluidProperties.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C
 $(forceModels)/Fines/FanningDynFines.C

src/recurrence/recModel/standardRecModel/standardRecModel.C

@@ -159,7 +159,24 @@ standardRecModel::standardRecModel
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 standardRecModel::~standardRecModel()
-{}
+{
+    const objectRegistry& objReg = base_.mesh().thisDb();
+
+    for(int j=0; j<volScalarFieldNames_.size(); j++)
+    {
+        objReg.checkOut(volScalarFieldList_[j][virtualTimeIndex]);
+    }
+
+    for(int j=0; j<volVectorFieldNames_.size(); j++)
+    {
+        objReg.checkOut(volVectorFieldList_[j][virtualTimeIndex]);
+    }
+
+    for(int j=0; j<surfaceScalarFieldNames_.size(); j++)
+    {
+        objReg.checkOut(surfaceScalarFieldList_[j][virtualTimeIndex]);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allclean.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+#------------------------------------------------------------------------------
+# Allclean script for falling sphere test case
+# clean CFD and DEM part
+# Achuth N. Balachandran Nair - October 2018
+#------------------------------------------------------------------------------
+
+#- source environment variables
+. ~/.bashrc
+
+#- source CFDEM functions
+source $CFDEM_PROJECT_DIR/etc/functions.sh
+
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+
+#- clean up case
+echo "deleting data at: $casePath"
+source $WM_PROJECT_DIR/bin/tools/CleanFunctions
+cd $casePath/CFD
+cleanCase
+rm -r $casePath/CFD/clockData
+rm $casePath/DEM/post/*.*
+#- preserve post directory
+touch $casePath/DEM/post/.gitignore
+echo "done"
+

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/Allrun.sh

@@ -0,0 +1,25 @@
+#!/bin/bash
+#------------------------------------------------------------------------------
+# allrun script for falling sphere testcase
+# run falling_sphere_two_way_coupling
+# Achuth N. Balachandran Nair - Oct. 2018
+#------------------------------------------------------------------------------
+
+source $CFDEM_PROJECT_DIR/etc/functions.sh
+
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+
+echo $casePath
+
+# 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 CFD/
+    blockMesh
+fi
+
+#- run parallel CFD-DEM
+bash $casePath/parCFDDEMrun.sh

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/U

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    channelWall
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/Us

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      Us;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    channelWall
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/p

@@ -1,38 +1,46 @@
 /*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      dSmoothing;
+    object      p;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 1 0 0 0 0 0];
+dimensions      [0 2 -2 0 0 0 0];
 
 internalField   uniform 0;
 
 boundaryField
 {
-    side-walls
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    channelWall
     {
         type            zeroGradient;
     }
 
-    inlet
-    {
-	type            zeroGradient;
-    }
-
     outlet
     {
-        type            zeroGradient;  
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    frontAndBack
+    {
+        type            empty;
     }
 }
+
 // ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/phiIB

@@ -1,17 +1,16 @@
 /*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
-    object      nut;
+    object      phiIB;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -21,24 +20,27 @@ internalField   uniform 0;
 
 boundaryField
 {
-    side-walls
+    inlet
     {
         type            fixedValue;
         value           uniform 0;
     }
 
-    inlet
+    channelWall
     {
-        type            calculated;
-        value           uniform 0;
+        type            zeroGradient;
     }
 
     outlet
     {
-        type            calculated;
+        type            fixedValue;
         value           uniform 0;
     }
+
+    frontAndBack
+    {
+        type            empty;
+    }
 }
 
-
 // ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/rho

@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
@@ -16,16 +16,16 @@ FoamFile
 
 dimensions      [1 -3 0 0 0 0 0];
 
-internalField   uniform 1.58;
+internalField   uniform 1;
 
 boundaryField
 {
-    side-walls
+    inlet
     {
         type            zeroGradient;
     }
 
-    inlet
+    channelWall
     {
         type            zeroGradient;
     }
@@ -36,5 +36,4 @@ boundaryField
     }
 }
 
-
 // ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/0/voidfraction

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      voidfraction;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    channelWall
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/couplingProperties

@@ -0,0 +1,101 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      couplingProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// sub-models & settings
+
+modelType none;
+
+checkPeriodicCells;
+
+verbose;
+
+couplingInterval 10;
+
+depth 0;
+
+voidFractionModel IB;
+
+locateModel engineIB;
+
+meshMotionModel noMeshMotion;
+
+regionModel allRegion;
+
+dataExchangeModel twoWayMPI;
+
+IOModel basicIO;
+
+probeModel off;
+
+averagingModel dilute;
+
+clockModel off;
+
+smoothingModel off;
+
+forceModels
+(
+    ShirgaonkarIB
+    ArchimedesIB
+);
+
+momCoupleModels
+(
+);
+
+turbulenceModelType "turbulenceProperties";
+
+// sub-model properties
+
+ShirgaonkarIBProps
+{
+    velFieldName "U";
+    pressureFieldName "p";
+    densityFieldName "rho";
+    verbose;
+    useTorque;
+}
+
+ArchimedesIBProps
+{
+    densityFieldName "rho";
+    gravityFieldName "g";
+    voidfractionFieldName "voidfractionNext";
+}
+
+twoWayMPIProps
+{
+    liggghtsPath "../DEM/in.liggghts_run";
+}
+
+IBProps
+{
+    maxCellsPerParticle 20000;
+    alphaMin 0.30;
+    scaleUpVol 1.0;
+}
+
+engineIBProps
+{
+    engineProps
+    {
+        treeSearch false;
+    }
+    zSplit 8;
+    xySplit 16;
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/dynamicMeshDict

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh       staticFvMesh;
+
+// ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/g

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+// ************************************************************************* //

tutorials/cfdemSolverIB/falling_sphere_two_way_coupling/CFD/constant/liggghtsCommands

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      liggghtsCommands;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+liggghtsCommandModels
+(
+   runLiggghts
+);
+
+// ************************************************************************* //