Merge pull request #115 from ParticulateFlow/release

Release 21.03

applications/solvers/cfdemSolverIB/Make/options

@@ -14,7 +14,8 @@ EXE_INC = \
     -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)/fvOptions/lnInclude \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/cfdemSolverMultiphase/Make/options

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

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/Make/options

@@ -8,7 +8,8 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -Wno-deprecated-copy
 
 LIB_LIBS = \
     -linterfaceProperties \

applications/solvers/cfdemSolverPiso/Make/options

@@ -10,6 +10,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/cfdemSolverPisoMS/Make/options

@@ -11,6 +11,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/cfdemSolverPisoScalar/Make/options

@@ -11,6 +11,7 @@ EXE_INC = \
     -I../cfdemSolverPiso \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/cfdemSolverRhoPimple/Make/options

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

applications/solvers/cfdemSolverRhoPimpleChem/Make/options

@@ -28,8 +28,7 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam \
-
-
+    -Wno-deprecated-copy
 
 
 EXE_LIBS = \

applications/solvers/rStatAnalysis/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/solvers/rcfdemSolverBase/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/solvers/rcfdemSolverCoupledHeattransfer/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/solvers/rcfdemSolverHeattransfer/Make/files

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

applications/solvers/rcfdemSolverHeattransfer/Make/options

@@ -1,27 +0,0 @@
-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$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
-    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
-    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-
-EXE_LIBS = \
-    -L$(CFDEM_LIB_DIR)\
-    -lrecurrence \
-    -lturbulenceModels \
-    -lincompressibleTurbulenceModels \
-    -lincompressibleTransportModels \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lfvOptions \
-    -l$(CFDEM_LIB_NAME) \
-     $(CFDEM_ADD_LIB_PATHS) \
-$(CFDEM_ADD_LIBS)

applications/solvers/rcfdemSolverHeattransfer/TEqImp.H

@@ -1,40 +0,0 @@
-    volScalarField rhoeps = rhoRec*voidfractionRec;
-
-    particleCloud.energyContributions(Qsource);
-
-    particleCloud.energyCoefficients(QCoeff);
-
-    K = 0.5*magSqr(URec);
-
-    addSource = fvc::ddt(rhoeps, K) + fvc::div(phiRec, K) +
-    fvc::div
-    (
-        fvc::absolute(phiRec/fvc::interpolate(rhoRec), voidfractionRec*URec), pRec
-    );
-
-
-    fvScalarMatrix TEqn =
-    (
-        fvm::ddt(rhoeps, T)
-        + fvm::div(phiRec, T)
-        + addSource/Cv
-        - fvm::laplacian(voidfractionRec*thCond/Cv, T)
-        - Qsource/Cv
-        - fvm::Sp(QCoeff/Cv, T)
-        ==
-        fvOptions(rhoeps, T)    // no fvOptions support yet
-    );
-
-    //TEqn.relax(relaxCoeff);
-
-    fvOptions.constrain(TEqn); // no fvOptions support yet
-
-    TEqn.solve();
-
-    particleCloud.clockM().start(31,"postFlow");
-    counter++;
-
-    if((counter - couplingSubStep) % dtDEM2dtCFD == 0)
-        particleCloud.postFlow();
-
-    particleCloud.clockM().stop("postFlow");

applications/solvers/rcfdemSolverHeattransfer/createFields.H

@@ -1,206 +0,0 @@
-    // dummy fields
-    Info << "\nCreating dummy pressure field\n" << endl;
-    volScalarField pRec
-    (
-        IOobject
-        (
-            "pRec",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    // recurrence fields
-    Info << "\nCreating recurrence fields.\n" << endl;
-
-    volScalarField rhoRec
-    (
-        IOobject
-        (
-            "rhoRec",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh//,
-        //dimensionedScalar("rhoRec", dimensionSet(1, -3, 0, 0, 0), 1.0)
-    );
-
-    volVectorField URec
-    (
-        IOobject
-        (
-            "URec",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField voidfractionRec
-    (
-        IOobject
-        (
-            "voidfractionRec",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volVectorField UsRec
-    (
-        IOobject
-        (
-            "UsRec",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    // heat transfer fields
-    Info << "\nCreating heat transfer fields.\n" << endl;
-
-    volScalarField Qsource
-    (
-        IOobject
-        (
-            "Qsource",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
-    );
-
-    volScalarField QCoeff
-    (
-        IOobject
-        (
-            "QCoeff",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
-    );
-
-    volScalarField thCond
-    (
-        IOobject
-        (
-            "thCond",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
-        "zeroGradient"
-    );
-
-    volScalarField T
-    (
-        IOobject
-        (
-            "T",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    // calculated fields
-    Info << "\nCreating fields subject to calculation\n" << endl;
-    volScalarField voidfraction
-    (
-        IOobject
-        (
-            "voidfraction",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        voidfractionRec
-    );
-
-    volVectorField Us
-    (
-        IOobject
-        (
-            "Us",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        UsRec
-    );
-
-    // write fields for t=t_start
-    voidfraction.write();
-    Us.write();
-//===============================
-
-
-    Info << "Calculating face flux field phiRec\n" << endl;
-    surfaceScalarField phiRec
-    (
-        IOobject
-        (
-            "phiRec",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        linearInterpolate(URec*voidfractionRec*rhoRec) & mesh.Sf()
-    );
-    phiRec.write();
-
-
-    singlePhaseTransportModel laminarTransport(URec, phiRec);
-
-    autoPtr<incompressible::turbulenceModel> turbulence
-    (
-        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
-    );
-
-    const IOdictionary& transportProps = mesh.lookupObject<IOdictionary>("transportProperties");
-    dimensionedScalar Cv(transportProps.lookup("Cv"));
-
-    volScalarField addSource
-    (
-        IOobject
-        (
-            "addSource",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
-    );
-
-    Info << "Creating field kinetic energy K\n" << endl;
-    volScalarField K("K", 0.5*magSqr(URec));

applications/solvers/rcfdemSolverHeattransfer/rcfdemSolverHeattransfer.C

@@ -1,128 +0,0 @@
-/*---------------------------------------------------------------------------*\
-    CFDEMcoupling academic - Open Source CFD-DEM coupling
-
-    Contributing authors:
-    Thomas Lichtenegger
-    Copyright (C) 2015- Johannes Kepler University, Linz
--------------------------------------------------------------------------------
-License
-    This file is part of CFDEMcoupling academic.
-
-    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    rcfdemSolverHeattransfer
-
-Description
-    Solves heat transfer between fluid and particles based on rCFD
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "singlePhaseTransportModel.H"
-#include "turbulentTransportModel.H"
-#include "fvOptions.H"
-
-#include "cfdemCloudRec.H"
-#include "recBase.H"
-#include "recModel.H"
-
-#include "cfdemCloudEnergy.H"
-#include "clockModel.H"
-#include "thermCondModel.H"
-#include "energyModel.H"
-
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createFields.H"
-    #include "createFvOptions.H"
-
-    cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
-    recBase recurrenceBase(mesh);
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
-
-    label recTimeIndex = 0;
-    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
-    scalar startTime = runTime.startTime().value();
-
-    // control coupling behavior in case of substepping
-    // assumes constant timestep size
-    label counter = 0;
-    label couplingSubStep = recurrenceBase.couplingSubStep();//5;//3;
-    double dtProp =  particleCloud.dataExchangeM().couplingTime() / runTime.deltaTValue();
-    label dtDEM2dtCFD = int(dtProp + 0.5);
-    Info << "deltaT_DEM / deltaT_CFD = " << dtDEM2dtCFD << endl;
-    if (dtDEM2dtCFD > 1)
-        Info << "coupling at substep " << couplingSubStep << endl;
-
-
-    while (runTime.run())
-    {
-        runTime++;
-
-        // do stuff (every lagrangian time step)
-        particleCloud.clockM().start(1,"Global");
-
-        Info << "Time = " << runTime.timeName() << nl << endl;
-
-        particleCloud.clockM().start(2,"Coupling");
-
-        particleCloud.evolve(voidfraction,Us,URec);
-
-        particleCloud.clockM().stop("Coupling");
-
-        particleCloud.clockM().start(26,"Flow");
-        #include "TEqImp.H"
-        particleCloud.clockM().stop("Flow");
-
-        particleCloud.clockM().start(32,"ReadFields");
-        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
-        {
-            recurrenceBase.updateRecFields();
-            #include "updateFields.H"
-            recTimeIndex++;
-        }
-        particleCloud.clockM().stop("ReadFields");
-
-        particleCloud.clockM().start(33,"Output");
-        runTime.write();
-        particleCloud.clockM().stop("Output");
-
-        particleCloud.clockM().stop("Global");
-
-        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-             << nl << endl;
-
-    }
-
-    Info << "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/rcfdemSolverHeattransfer/updateFields.H

@@ -1,14 +0,0 @@
-recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
-recurrenceBase.recM().exportVolScalarField("rho",rhoRec);
-recurrenceBase.recM().exportVolScalarField("p",pRec);
-recurrenceBase.recM().exportVolVectorField("Us",UsRec);
-recurrenceBase.recM().exportVolVectorField("U",URec);
-recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);
-
-{
-    volScalarField& NuField(const_cast<volScalarField&>(mesh.lookupObject<volScalarField> ("NuField")));
-    recurrenceBase.recM().exportVolScalarField("NuField",NuField);
-}
-
-
-#include "updateRho.H"

applications/solvers/rcfdemSolverHeattransfer/updateRho.H

@@ -1,32 +0,0 @@
-// work-around for transient properties
-// needs to be specified for each case
-
-// case 1
-
-forAll(rhoRec,cellI)
-{
-    if (mesh.C()[cellI].z() < 0.00228)
-        rhoRec[cellI] = 1.18+(1.085-1.18)*Foam::exp(-0.065*runTime.timeOutputValue());
-    else if (mesh.C()[cellI].z() < 0.00456)
-        rhoRec[cellI] = 1.18+(1.01-1.18)*Foam::exp(-0.05*runTime.timeOutputValue());
-    else if (mesh.C()[cellI].z() < 0.00684)
-        rhoRec[cellI] = 1.18+(0.98-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());
-    else
-        rhoRec[cellI] = 1.18+(0.955-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());   
-}
-
-
-// case 2
-/*
-forAll(rhoRec,cellI)
-{
-    if (mesh.C()[cellI].z() < 0.00228)
-        rhoRec[cellI] = 1.18+(1.115-1.18)*Foam::exp(-0.065*runTime.timeOutputValue());
-    else if (mesh.C()[cellI].z() < 0.00456)
-        rhoRec[cellI] = 1.18+(1.04-1.18)*Foam::exp(-0.05*runTime.timeOutputValue());
-    else if (mesh.C()[cellI].z() < 0.00684)
-        rhoRec[cellI] = 1.18+(1.005-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());
-    else
-        rhoRec[cellI] = 1.18+(0.96-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());   
-}
-*/

applications/solvers/rcfdemSolverRhoSteadyPimple/Make/options

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

applications/solvers/rtfmSolverSpecies/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/cfdemPostproc/Make/options

@@ -9,6 +9,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -Wno-deprecated-copy
 
 
 EXE_LIBS = \

applications/utilities/rBaseMirror/rBaseMirrorScalar/Make/options

@@ -7,7 +7,8 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/fvOptions/lnInclude 
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/options

@@ -7,7 +7,8 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/fvOptions/lnInclude 
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/utilities/rSmoother/Make/options

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

doc/CFDEMcoupling_solvers.txt

@@ -18,7 +18,6 @@ This section lists all CFDEMcoupling solvers alphabetically.
 rStatAnalysis,
 rcfdemSolverBase,
 rcfdemSolverCoupledHeattransfer,
-rcfdemSolverHeattransfer,
 rcfdemSolverRhoSteadyPimple,
 recSolverTurbTransport,
 rtfmSolverSpecies,

doc/forceModel_particleCellVolume.txt

@@ -22,11 +22,14 @@ particleCellVolumeProps
     lowerThreshold scalar2;
     startTime      scalar3;
     verbose;
+    writeToFile    switch1;
 \} :pre
 
 {scalar1} = only cells with a field value (magnitude) lower than this upper threshold are considered :l
 {scalar2} = only cells with a field value (magnitude) greater than this lower threshold are considered :l
 {scalar3} = (optional, default 0) start time of volume calculation and output :l
+{verbose} = (optional, default false) keyword only (mostly used for debugging) :l
+{switch1} = (optional, default false) switch for file output :l
 :ule
 
 [Examples:]

doc/forceModel_volWeightedAverage.txt

@@ -30,6 +30,7 @@ volWeightedAverageProps
     upperThreshold     scalar1;
     lowerThreshold     scalar2;
     verbose;
+    writeToFile        switch1;
 \} :pre
 
 {time} = (optional, default 0.) time to start the averaging :ulb,l
@@ -38,6 +39,7 @@ volWeightedAverageProps
 {scalar1} = only cells with a field value (magnitude) lower than this upper threshold are considered :l
 {scalar2} = only cells with a field value (magnitude) greater than this lower threshold are considered :l
 {verbose} = (optional, default false) keyword only (mostly used for debugging) :l
+{switch1} = (optional, default false) switch for file output :l
 :ule
 
 [Examples:]

etc/bashrc

@@ -17,7 +17,7 @@
 #------------------------------------------------------------------------------
 
 export CFDEM_PROJECT=CFDEM
-export CFDEM_VERSION=20.05
+export CFDEM_VERSION=21.03
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
@@ -199,6 +199,9 @@ alias cfdemCleanCFDEM='bash $CFDEM_PROJECT_DIR/etc/cleanCFDEMcoupling.sh'
 #- shortcut to compile LIGGGHTS + sublibraries
 alias cfdemCompLIG='bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS.sh'
 
+#- shortcut to compile LIGGGHTS sublibraries
+alias cfdemCompLIGlib='bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS_lib.sh'
+
 #- shortcut to compile CFDEMcoupling +LIGGGHTS
 alias cfdemCompCFDEMall='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_all.sh'
 

etc/cshrc

@@ -15,7 +15,7 @@
 #------------------------------------------------------------------------------
 
 setenv CFDEM_PROJECT CFDEM
-setenv CFDEM_VERSION 20.05
+setenv CFDEM_VERSION 21.03
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
@@ -232,6 +232,9 @@ alias cfdemCleanCFDEM 'bash $CFDEM_PROJECT_DIR/etc/cleanCFDEMcoupling.sh'
 #- shortcut to compile LIGGGHTS + sublibraries
 alias cfdemCompLIG 'bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS.sh'
 
+#- shortcut to compile LIGGGHTS sublibraries
+alias cfdemCompLIGlib 'bash $CFDEM_PROJECT_DIR/etc/compileLIGGGHTS_lib.sh'
+
 #- shortcut to compile CFDEMcoupling +LIGGGHTS
 alias cfdemCompCFDEMall 'bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_all.sh'
 

etc/solver-list.txt

@@ -1,5 +1,4 @@
 rcfdemSolverRhoSteadyPimple/dir
-rcfdemSolverHeattransfer/dir
 rcfdemSolverCoupledHeattransfer/dir
 rStatAnalysis/dir
 rcfdemSolverBase/dir

src/finiteVolume/Make/options

@@ -4,7 +4,8 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -Wno-old-style-cast \
     -Wno-unused-result \
-    -Wno-literal-suffix
+    -Wno-literal-suffix \
+    -Wno-deprecated-copy
 
 LIB_LIBS = \
     -lOpenFOAM \

src/lagrangian/cfdemParticle/Make/files

@@ -91,6 +91,7 @@ $(forceModels)/pdCorrelation/pdCorrelation.C
 $(forceModels)/surfaceTensionForce/surfaceTensionForce.C
 $(forceModels)/gradPForceSmooth/gradPForceSmooth.C
 $(forceModels)/particleDeformation/particleDeformation.C
+$(forceModels)/turbulentDispersion/turbulentDispersion.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C

src/lagrangian/cfdemParticle/Make/options

@@ -23,7 +23,8 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/cfdTools \
     -Wno-old-style-cast \
     -Wno-unused-result \
-    -Wno-literal-suffix
+    -Wno-literal-suffix \
+    -Wno-deprecated-copy
 
 
 LIB_LIBS = \

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

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

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -84,10 +84,12 @@ cfdemCloud::cfdemCloud
     ignore_(couplingProperties_.found("ignore")),
     allowCFDsubTimestep_(true),
     limitDEMForces_(couplingProperties_.found("limitDEMForces")),
+    phaseInForces_(couplingProperties_.found("phaseInForcesTime")),
     getParticleDensities_(couplingProperties_.lookupOrDefault<bool>("getParticleDensities",false)),
     getParticleEffVolFactors_(couplingProperties_.lookupOrDefault<bool>("getParticleEffVolFactors",false)),
     getParticleTypes_(couplingProperties_.lookupOrDefault<bool>("getParticleTypes",false)),
     maxDEMForce_(0.),
+    phaseInForcesTime_(couplingProperties_.lookupOrDefault<scalar>("phaseInForcesTime",0.0)),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
     velocities_(NULL),
@@ -158,6 +160,7 @@ cfdemCloud::cfdemCloud
             turbulenceModelType_
         )
     ),
+    particlePropertyTable(32),
     dataExchangeModel_
     (
         dataExchangeModel::New
@@ -368,6 +371,12 @@ cfdemCloud::cfdemCloud
         if (verbose_) Info << "nPatchesNonCyclic=" << nPatchesNonCyclic << ", nPatchesCyclic=" << nPatchesCyclic << endl;
         Warning << "Periodic handing is disabled because the domain is not fully periodic!\n" << endl;
     }
+
+    // check if phasing-in time existing and is meaningful
+    if (phaseInForces_ && phaseInForcesTime_ < SMALL)
+    {
+        FatalError << "phasing-in time too small" << endl;
+    }
 }
 
 // * * * * * * * * * * * * * * * * Destructors  * * * * * * * * * * * * * * //
@@ -394,6 +403,22 @@ cfdemCloud::~cfdemCloud()
     if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1);
     if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1);
     if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1);
+
+    for
+    (
+        HashTable<particleProperty>::iterator iter = particlePropertyTable.begin();
+        iter != particlePropertyTable.end();
+        ++iter
+    )
+    {
+        if ((*(iter().ti)) == typeid(int**)) {
+            dataExchangeM().destroy(iter().ref<int**>(),-1);
+        } else if ((*(iter().ti)) == typeid(double**)) {
+            dataExchangeM().destroy(iter().ref<double**>(),-1);
+        } else {
+            FatalError << "Trying to destroy property of type " << iter().ti->name() << endl;
+        }
+    }
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -468,6 +493,19 @@ void cfdemCloud::setForces()
         }
         Info << "largest particle-fluid interaction on particle: " << maxF << endl;
     }
+
+    if (phaseInForces_)
+    {
+        scalar tfrac = (mesh_.time().timeOutputValue()-mesh_.time().startTime().value())/phaseInForcesTime_;
+        if (tfrac <= 1.0)
+        {
+            for (int index = 0;index <  numberOfParticles(); ++index)
+            {
+                for(int i=0;i<3;i++) DEMForces_[index][i] *= tfrac;
+                Cds_[index][0] *= tfrac;
+            }
+        }
+    }
 }
 
 void cfdemCloud::setParticleForceField()
@@ -759,9 +797,60 @@ bool cfdemCloud::reAllocArrays()
         if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1);
         if(getParticleEffVolFactors_) dataExchangeM().allocateArray(particleEffVolFactors_,0.,1);
         if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0,1);
+
+        for
+        (
+            HashTable<particleProperty>::iterator iter = particlePropertyTable.begin();
+            iter != particlePropertyTable.end();
+            ++iter
+        )
+        {
+            if (iter().size > 0) {
+                ///Info << "!! about to realloc property of type " << iter().ti->name() << endl;
+                if ((*(iter().ti)) == typeid(int**)) {
+                    dataExchangeM().allocateArray(iter().ref<int**>(),iter().initVal,iter().size);
+                } else if ((*(iter().ti)) == typeid(double**)) {
+                    dataExchangeM().allocateArray(iter().ref<double**>(),iter().initVal,iter().size);
+                } else {
+                    FatalError << "Trying to realloc property of type " << iter().ti->name() << endl;
+                }
+            }
+        }
+
         arraysReallocated_ = true;
         return true;
     }
+    else
+    {
+        for
+        (
+            HashTable<particleProperty>::iterator iter = particlePropertyTable.begin();
+            iter != particlePropertyTable.end();
+            ++iter
+        )
+        {
+            if (iter().size > 0 && iter().reset) {
+                if ((*(iter().ti)) == typeid(int**)) {
+                    int**& property = iter().ref<int**>();
+                    for (int index=0; index<numberOfParticles(); ++index) {
+                        for (int icomponent=0; icomponent<iter().size; ++icomponent) {
+                            property[index][icomponent] = iter().initVal;
+                        }
+                    }
+                } else if ((*(iter().ti)) == typeid(double**)) {
+                    double**& property = iter().ref<double**>();
+                    for (int index=0; index<numberOfParticles(); ++index) {
+                        for (int icomponent=0; icomponent<iter().size; ++icomponent) {
+                            property[index][icomponent] = iter().initVal;
+                        }
+                    }
+                } else {
+                    FatalError << "Trying to reset property of type " << iter().ti->name() << endl;
+                }
+            }
+        }
+    }
+
     return false;
 }
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -45,6 +45,7 @@ SourceFiles
 // choose version
 #include "OFversion.H"
 #include <vector>
+#include <typeinfo>
 
 #include "fvCFD.H"
 #include "IFstream.H"
@@ -96,6 +97,8 @@ protected:
 
     const bool limitDEMForces_;
 
+    const bool phaseInForces_;
+
     const bool getParticleDensities_;
 
     const bool getParticleEffVolFactors_;
@@ -104,6 +107,8 @@ protected:
 
     scalar maxDEMForce_;
 
+    scalar phaseInForcesTime_;
+
     const word modelType_;
 
     double **positions_; // particle positions
@@ -184,6 +189,21 @@ protected:
 
     const turbulenceModel& turbulence_;
 
+    struct particleProperty {
+        void** property;            // pointer to per-particle data; memory deallocation by cfdemCloud
+        const std::type_info* ti;   // type of particle data (int**, double**)
+        int size;                   // size of single particle data; memory allocation by cfdemCloud if size > 0
+        double initVal;             // initial property value
+        bool reset;                 // if true, data is reset to initial value every coupling step
+        template<typename T>
+        T& ref() {
+          if (*ti == typeid(T)) return *reinterpret_cast<T*>(&property);
+          else throw std::bad_cast();
+        }
+    };
+
+    HashTable<particleProperty> particlePropertyTable; // table of registered per-particle properties
+
     autoPtr<dataExchangeModel> dataExchangeModel_;
 
     PtrList<forceModel> forceModel_;
@@ -439,6 +459,14 @@ public:
         void otherForces(volVectorField&);
 
         bool checkPeriodicCells() const { return checkPeriodicCells_; }
+
+        template<typename T>
+        void registerParticleProperty(const word& property, int size=0, double initVal=0.0, bool reset=true);
+        template<typename T>
+        T& getParticlePropertyRef(const word& property);
+protected:
+        virtual int**& getParticlePropertyImpl(const word& property, int**);
+        virtual double**& getParticlePropertyImpl(const word& property, double**);
 };
 
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -403,5 +403,28 @@ inline const turbulenceModel& cfdemCloud::turbulence() const
     return turbulence_;
 }
 
+
+template<typename T>
+void cfdemCloud::registerParticleProperty(const word& property, int size, double initVal, bool reset)
+{
+    particlePropertyTable.insert(property,{NULL,&typeid(T),size,initVal,reset});
+}
+
+template<typename T>
+T& cfdemCloud::getParticlePropertyRef(const word& property)
+{
+    return getParticlePropertyImpl(property, static_cast<T>(0));
+}
+
+inline int**& cfdemCloud::getParticlePropertyImpl(const word& property, int**)
+{
+    return particlePropertyTable[property].ref<int**>();
+}
+
+inline double**& cfdemCloud::getParticlePropertyImpl(const word& property, double**)
+{
+    return particlePropertyTable[property].ref<double**>();
+}
+
 }
 // ************************************************************************* //

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

@@ -71,15 +71,17 @@ diffusionCoefficient::diffusionCoefficient
     pressureFieldName_(propsDict_.lookupOrDefault<word>("pressureFieldName","p")),
     P_(sm.mesh().lookupObject<volScalarField>(pressureFieldName_)),
     partPressureName_(propsDict_.lookupOrDefault<word>("partPressureName","partP")),
-    partPressure_(NULL),
     X_(speciesNames_.size()),
     diffusantGasNames_(propsDict_.lookup("diffusantGasNames")),
-    diffusionCoefficients_(diffusantGasNames_.size(),NULL),
     Xdiffusant_(diffusantGasNames_.size()),
     initialized_(false)
 {
     particleCloud_.checkCG(false);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partPressureName_,1);
+    for (int i=0; i<diffusantGasNames_.size(); i++)
+    {
+        particleCloud_.registerParticleProperty<double**>(diffusantGasNames_[i],1);
+    }
     createCoeffs();
     molWeightTable();
 }
@@ -88,37 +90,10 @@ diffusionCoefficient::diffusionCoefficient
 
 diffusionCoefficient::~diffusionCoefficient()
 {
-    particleCloud_.dataExchangeM().destroy(partPressure_,1);
-    for (int i=0; i<diffusantGasNames_.size(); i++) particleCloud_.dataExchangeM().destroy(diffusionCoefficients_[i],1);
-
-    coeffs.clearStorage();
-    molWeight.clearStorage();
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-void diffusionCoefficient::allocateMyArrays() const
-{
-    double initVal=0.0;
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partPressure_,initVal,1,"nparticles");
-        for (int i=0; i<diffusantGasNames_.size(); i++)
-        {
-            particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1,"nparticles");
-        }
-    }
-}
-
-void diffusionCoefficient::reAllocMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partPressure_,initVal,1,"nparticles");
-    for (int i=0; i<diffusantGasNames_.size(); i++)
-    {
-        particleCloud_.dataExchangeM().allocateArray(diffusionCoefficients_[i],initVal,1);
-    }
-}
 
 void diffusionCoefficient::init()
 {
@@ -156,9 +131,6 @@ void diffusionCoefficient::execute()
         init();
     }
 
-    // realloc the arrays
-    reAllocMyArrays();
-
     label  cellI=0;
     scalar Tfluid(0);
     scalar Pfluid(0);
@@ -172,6 +144,8 @@ void diffusionCoefficient::execute()
     interpolationCellPoint <scalar> TInterpolator_(tempField_);
     interpolationCellPoint <scalar> PInterpolator_(P_);
 
+    double**& partPressure_ = particleCloud_.getParticlePropertyRef<double**>(partPressureName_);
+
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
@@ -207,6 +181,7 @@ void diffusionCoefficient::execute()
 
             for (int j=0; j<diffusantGasNames_.size(); j++)
             {
+                double**& diffusionCoefficients_ = particleCloud_.getParticlePropertyRef<double**>(diffusantGasNames_[j]);
                 TotalFraction_[j] = 0.0;
                 dBinary_ = 0.0;
 
@@ -243,9 +218,9 @@ void diffusionCoefficient::execute()
 
                             // pass on dCoeff values to array
                             if (TotalFraction_[j] < VSMALL)
-                                diffusionCoefficients_[j][index][0] = VSMALL;
+                                diffusionCoefficients_[index][0] = VSMALL;
                             else
-                                diffusionCoefficients_[j][index][0] = (1.0 - Xdiffusant_[j][cellI]) / TotalFraction_[j];
+                                diffusionCoefficients_[index][0] = (1.0 - Xdiffusant_[j][cellI]) / TotalFraction_[j];
                         }
                         else
                         {
@@ -258,7 +233,7 @@ void diffusionCoefficient::execute()
                 }
 
                 if(verbose_)
-                    Info << "diffusionCoefficient of species " << diffusantGasNames_[j] << " = " << diffusionCoefficients_[j][index][0] << endl;
+                    Info << "diffusionCoefficient of species " << diffusantGasNames_[j] << " = " << diffusionCoefficients_[index][0] << endl;
             }
         }
     }
@@ -268,7 +243,8 @@ void diffusionCoefficient::execute()
     for (int j=0; j<diffusantGasNames_.size(); j++)
     {
         word pushName = diffusantGasNames_[j] + "_diffCoeff";
-        particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_[j]);
+        double**& diffusionCoefficients_ = particleCloud_.getParticlePropertyRef<double**>(diffusantGasNames_[j]);
+        particleCloud_.dataExchangeM().giveData(pushName,"scalar-atom",diffusionCoefficients_);
     }
 
     Info << "give data done" << endl;

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

@@ -76,19 +76,15 @@ private:
 
     word partPressureName_;
 
-    mutable double **partPressure_;
-
     UPtrList<volScalarField> X_;
 
     wordList diffusantGasNames_;
 
-    mutable List<double**> diffusionCoefficients_;
-
     UPtrList<volScalarField> Xdiffusant_;
 
-    HashTable<scalar, word> coeffs;
+    HashTable<scalar> coeffs;
 
-    HashTable<scalar, word> molWeight;
+    HashTable<scalar> molWeight;
 
     void createCoeffs();
 
@@ -100,8 +96,6 @@ private:
     // calculate denominator part diffusion volume equation
     double calcDiffVol(int, int);
 
-    void allocateMyArrays() const;
-
     bool initialized_;
 
     void init();
@@ -126,8 +120,6 @@ private:
 
     // Member Functions
         void execute();
-
-        void reAllocMyArrays() const;
 };
 
 

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

@@ -63,51 +63,26 @@ massTransferCoeff::massTransferCoeff
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partNuName_(propsDict_.lookupOrDefault<word>("partViscos","partNu")),
-    partNu_(NULL),
-    partReynolds_(propsDict_.lookupOrDefault<word>("partReynolds","partRe")),
-    partRe_(NULL),
+    partReName_(propsDict_.lookupOrDefault<word>("partReynolds","partRe")),
     scaleDia_(1)
 {
     particleCloud_.checkCG(true);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partNuName_,1);
+    particleCloud_.registerParticleProperty<double**>(partReName_,1);
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 massTransferCoeff::~massTransferCoeff()
 {
-    int nP_ = particleCloud_.numberOfParticles();
-
-    particleCloud_.dataExchangeM().destroy(partNu_,nP_);
-    particleCloud_.dataExchangeM().destroy(partRe_,nP_);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void massTransferCoeff::allocateMyArrays() const
-{
-    double initVal=0.0;
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d arrays
-        particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1,"nparticles");
-        particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1,"nparticles");
-    }
-}
-
-void massTransferCoeff::reAllocMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
-    particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void massTransferCoeff::execute()
 {
-    // realloc the arrays
-    reAllocMyArrays();
-
     #ifdef compre
         const volScalarField nufField = particleCloud_.turbulence().mu()/rho_;
     #else
@@ -139,6 +114,9 @@ void massTransferCoeff::execute()
     interpolationCellPoint <vector> UluidInterpolator_(U_);
     interpolationCellPoint <scalar> voidfractionInterpolator_(voidfraction_);
 
+    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuName_);
+    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReName_);
+
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
         cellI=particleCloud_.cellIDs()[index][0];
@@ -195,7 +173,7 @@ void massTransferCoeff::execute()
 
     // give DEM data
     particleCloud_.dataExchangeM().giveData(partNuName_, "scalar-atom", partNu_);
-    particleCloud_.dataExchangeM().giveData(partReynolds_, "scalar-atom", partRe_);
+    particleCloud_.dataExchangeM().giveData(partReName_, "scalar-atom", partRe_);
 
     Info << "give data done" << endl;
 }

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

@@ -56,27 +56,21 @@ private:
 
     const fvMesh& mesh_;
 
-    word velFieldName_;
+    const word velFieldName_;
 
     const volVectorField& U_;
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
-    word densityFieldName_;
+    const word densityFieldName_;
 
     const volScalarField& rho_;
 
-    word partNuName_;
+    const word partNuName_;
 
-    mutable double **partNu_;
-
-    word partReynolds_;
-
-    mutable double **partRe_;
-
-    void allocateMyArrays() const;
+    const word partReName_;
 
     mutable scalar scaleDia_;
 
@@ -101,8 +95,6 @@ public:
     // Member Functions
         void execute();
 
-        void reAllocMyArrays() const;
-
 };
 
 

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

@@ -57,7 +57,7 @@ reactantPerParticle::reactantPerParticle
     propsDict_(dict.subDict(typeName + "Props")),
     mesh_(sm.mesh()),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
-    reactantPerParticle_(NULL),
+    partReactantName_("reactantPerParticle"),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_)),
     particlesPerCell_
@@ -76,32 +76,16 @@ reactantPerParticle::reactantPerParticle
     Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1))
 {
     particleCloud_.checkCG(false);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partReactantName_,1);
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 reactantPerParticle::~reactantPerParticle()
 {
-    particleCloud_.dataExchangeM().destroy(reactantPerParticle_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void reactantPerParticle::allocateMyArrays() const
-{
-    double initVal=0.0;
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d arrays
-        particleCloud_.dataExchangeM().allocateArray(reactantPerParticle_,initVal,1,"nparticles");
-    }
-}
-
-void reactantPerParticle::reAllocMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(reactantPerParticle_,initVal,1);
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
@@ -112,8 +96,6 @@ void reactantPerParticle::execute()
     {
         return;
     }
-    // realloc the arrays
-    reAllocMyArrays();
 
     particlesPerCell_ *= 0.0;
 
@@ -121,6 +103,7 @@ void reactantPerParticle::execute()
     scalar voidfraction(1);
     scalar cellvolume(0.0);
     scalar particlesPerCell(1.0);
+    double**& reactantPerParticle_ = particleCloud_.getParticlePropertyRef<double**>(partReactantName_);
 
     // first create particles per cell field
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
@@ -147,10 +130,10 @@ void reactantPerParticle::execute()
         if (verbose_) Info << "reactantPerParticle_" << reactantPerParticle_[index][0] << endl;
     }
 
-        // give DEM data
-        particleCloud_.dataExchangeM().giveData("reactantPerParticle", "scalar-atom", reactantPerParticle_);
+    // give DEM data
+    particleCloud_.dataExchangeM().giveData(partReactantName_, "scalar-atom", reactantPerParticle_);
 
-        Info << "give data done" << endl;
+    Info << "give data done" << endl;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -56,16 +56,14 @@ private:
 
     bool verbose_;
 
-    mutable double **reactantPerParticle_;
+    const word partReactantName_;
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
     mutable volScalarField particlesPerCell_;
 
-    void allocateMyArrays() const;
-
     label loopCounter_;
 
     label Nevery_;
@@ -91,8 +89,6 @@ public:
     // Member Functions
         void execute();
 
-        void reAllocMyArrays() const;
-
 };
 
 

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

@@ -70,8 +70,6 @@ species::species
     speciesNames_(specDict_.lookup("species")),
     mod_spec_names_(speciesNames_.size()),
     X_(speciesNames_.size()),                           //volumeScalarFields of molarFractions
-    molarFractions_(speciesNames_.size(),NULL),         //the value of molar fractions for every species
-    changeOfSpeciesMass_(speciesNames_.size(),NULL),    //the values that are received from DEM with the name of Modified_+species name
     changeOfSpeciesMassFields_(speciesNames_.size()),   //the scalar fields generated with the values from Modified_+species names
     changeOfGasMassField_                               //the total change of Gas Mass field (when the Modified species
     (
@@ -89,18 +87,15 @@ species::species
     tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
     tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
     partTempName_(propsDict_.lookupOrDefault<word>("partTempName","partTemp")),
-    partTemp_(NULL),
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partRhoName_(propsDict_.lookupOrDefault<word>("partRhoName","partRho")),
-    partRho_(NULL),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_)),
     // total mole field
     molarConcFieldName_(propsDict_.lookupOrDefault<word>("totalMoleFieldName","molarConc")),
     molarConc_(sm.mesh().lookupObject<volScalarField>(molarConcFieldName_)),
     partMolarConcName_(propsDict_.lookupOrDefault<word>("partMoleName","partMolarConc")),
-    partMolarConc_(NULL),
     loopCounter_(-1),
     Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1)),
     massSourceCurr_(0.0),
@@ -108,56 +103,25 @@ species::species
     initialized_(false)
 {
     particleCloud_.checkCG(false);
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_,1);
+    particleCloud_.registerParticleProperty<double**>(partRhoName_,1);
+    particleCloud_.registerParticleProperty<double**>(partMolarConcName_,1);
+
+    for (int i=0; i<speciesNames_.size(); i++)
+    {
+        particleCloud_.registerParticleProperty<double**>("X_"+speciesNames_[i],1);
+        particleCloud_.registerParticleProperty<double**>("Modified_"+speciesNames_[i],1);
+    }
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 species::~species()
 {
-    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    particleCloud_.dataExchangeM().destroy(partRho_,1);
-    particleCloud_.dataExchangeM().destroy(partMolarConc_,1);
-
-
-    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(molarFractions_[i],1);
-    for (int i=0; i<speciesNames_.size();i++) particleCloud_.dataExchangeM().destroy(changeOfSpeciesMass_[i],1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-void species::allocateMyArrays() const
-{
-    double initVal=0.0;
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d arrays
-        particleCloud_.dataExchangeM().allocateArray(partRho_,initVal,1,"nparticles");
-        particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1,"nparticles");
-        particleCloud_.dataExchangeM().allocateArray(partMolarConc_,initVal,1,"nparticles");
-
-        for (int i=0; i<speciesNames_.size(); i++)
-        {
-            particleCloud_.dataExchangeM().allocateArray(molarFractions_[i],initVal,1,"nparticles");
-            particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1,"nparticles");
-        }
-    }
-}
-
-void species::reAllocMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partRho_,initVal,1);
-    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);
-    particleCloud_.dataExchangeM().allocateArray(partMolarConc_,initVal,1);
-
-    for (int i=0; i<speciesNames_.size(); i++)
-    {
-        particleCloud_.dataExchangeM().allocateArray(molarFractions_[i],initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1);
-    }
-}
-
 void species::init()
 {
     if(verbose_)
@@ -219,8 +183,6 @@ void species::execute()
     {
         return;
     }
-    // realloc the arrays
-    reAllocMyArrays();
 
     // get X_i, T, rho at particle positions
     label  cellI = 0;
@@ -235,6 +197,9 @@ void species::execute()
     interpolationCellPoint <scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint <scalar> molarConcInterpolator_(molarConc_);
 
+    double**& partRho_ = particleCloud_.getParticlePropertyRef<double**>(partRhoName_);
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partMolarConc_ = particleCloud_.getParticlePropertyRef<double**>(partMolarConcName_);
 
     for (int index=0; index<particleCloud_.numberOfParticles(); ++index)
     {
@@ -264,7 +229,8 @@ void species::execute()
             for (int i=0; i<speciesNames_.size();i++)
             {
                 // attention for indices when not communicating all species
-                molarFractions_[i][index][0] = X_[i][cellI];
+                double**& molarFractions_ = particleCloud_.getParticlePropertyRef<double**>("X_"+speciesNames_[i]);
+                molarFractions_[index][0] = X_[i][cellI];
             }
         }
     }
@@ -274,8 +240,9 @@ void species::execute()
     {
         for(int i =0; i<speciesNames_.size(); i++)
         {
+            double**& molarFractions_ = particleCloud_.getParticlePropertyRef<double**>("X_"+speciesNames_[i]);
             Info << "X_i = " << X_[i].name() << endl;
-            Info << "molarFractions_= " << molarFractions_[i][0][0] << endl;
+            Info << "molarFractions_= " << molarFractions_[0][0] << endl;
             Info << "partRho_[index][0] = " << partRho_[0][0] << endl;
             Info << "rhofluid = " << rhofluid << endl;
             Info << "partTemp_[index][0] = " << partTemp_[0][0] << endl;
@@ -292,7 +259,8 @@ void species::execute()
 
         for (int i=0; i<speciesNames_.size();i++)
         {
-            particleCloud_.dataExchangeM().giveData("X_"+speciesNames_[i],"scalar-atom",molarFractions_[i]);
+            double**& molarFractions_ = particleCloud_.getParticlePropertyRef<double**>("X_"+speciesNames_[i]);
+            particleCloud_.dataExchangeM().giveData("X_"+speciesNames_[i],"scalar-atom",molarFractions_);
         }
 
         if (verbose_) Info << "give data done" << endl;
@@ -305,17 +273,18 @@ void species::execute()
         changeOfGasMassField_.boundaryFieldRef() = 0.0;
         for (int i=0; i<speciesNames_.size();i++)
         {
+            double**& changeOfSpeciesMass_ = particleCloud_.getParticlePropertyRef<double**>("Modified_"+speciesNames_[i]);
             changeOfSpeciesMassFields_[i].primitiveFieldRef() = 0.0;
             changeOfSpeciesMassFields_[i].boundaryFieldRef() = 0.0;
 
-            particleCloud_.dataExchangeM().getData(mod_spec_names_[i],"scalar-atom",changeOfSpeciesMass_[i],particleCloud_.dataExchangeM().couplingInterval());
+            particleCloud_.dataExchangeM().getData(mod_spec_names_[i],"scalar-atom",changeOfSpeciesMass_,particleCloud_.dataExchangeM().couplingInterval());
 
-            if (verbose_) Info << "changeOfSpeciesMass received from DEM = " << changeOfSpeciesMass_[i][0][0] << endl;
+            if (verbose_) Info << "changeOfSpeciesMass received from DEM = " << changeOfSpeciesMass_[0][0] << endl;
 
             particleCloud_.averagingM().setScalarSumCentre
             (
                 changeOfSpeciesMassFields_[i],
-                changeOfSpeciesMass_[i],
+                changeOfSpeciesMass_,
                 particleCloud_.particleWeights(),
                 NULL
             );

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

@@ -68,44 +68,32 @@ private:
 
     UPtrList<volScalarField> X_;
 
-    mutable List<double**> molarFractions_;
-
-    mutable List<double**> changeOfSpeciesMass_;
-
     PtrList<volScalarField> changeOfSpeciesMassFields_;
 
     volScalarField changeOfGasMassField_;
 
-    word tempFieldName_;
+    const word tempFieldName_;
 
     const volScalarField& tempField_;              // ref to gas temperature field
 
-    word partTempName_;
+    const word partTempName_;
 
-    mutable double **partTemp_;                    // gas temperature at particle positions
-
-    word densityFieldName_;
+    const word densityFieldName_;
 
     const volScalarField& rho_;
 
-    word partRhoName_;
+    const word partRhoName_;
 
-    mutable double **partRho_;                     // gas density at particle positions
-
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
-    void allocateMyArrays() const;
-
     // total mole field
-    word molarConcFieldName_;
+    const word molarConcFieldName_;
 
     const volScalarField& molarConc_;
 
-    word partMolarConcName_;
-
-    mutable double **partMolarConc_;
+    const word partMolarConcName_;
 
     label loopCounter_;
 
@@ -141,8 +129,6 @@ public:
     // Member Functions
         void execute();
 
-        void reAllocMyArrays() const;
-
         tmp <volScalarField> Smi(const label i) const;
 
         tmp <volScalarField> Sm() const;

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

@@ -129,31 +129,31 @@ public:
         template <typename T>
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             T ** const& field
         ) const { getData(name,type,field,couplingStep_-1); }
 
         virtual void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const = 0;
 
         virtual void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const = 0;
 
         virtual void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype="double"
         ) const = 0;

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

@@ -78,24 +78,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const {}
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const {}
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype = ""
         ) const {}

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

@@ -84,8 +84,8 @@ oneWayVTK::~oneWayVTK()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 void oneWayVTK::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label step
 ) const
@@ -194,8 +194,8 @@ void oneWayVTK::getData
 
 void oneWayVTK::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* datatype
 ) const

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

@@ -87,16 +87,16 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const
@@ -104,8 +104,8 @@ public:
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype = ""
         ) const;

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

@@ -80,7 +80,7 @@ twoWayFiles::~twoWayFiles()
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-fileName twoWayFiles::getFilePath(word& name, bool in) const
+fileName twoWayFiles::getFilePath(const word& name, bool in) const
 {
     char timeStep[40];
 
@@ -100,7 +100,7 @@ fileName twoWayFiles::getFilePath(word& name, bool in) const
     return particleFilePathOld;
 }
 
-void twoWayFiles::renameFilePath(fileName& particleFilePathOld,word& name) const
+void twoWayFiles::renameFilePath(const fileName& particleFilePathOld, const word& name) const
 {
     char timeStep[40];
 
@@ -116,8 +116,8 @@ void twoWayFiles::renameFilePath(fileName& particleFilePathOld,word& name) const
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 void twoWayFiles::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label step
 ) const
@@ -164,8 +164,8 @@ void twoWayFiles::getData
 
 void twoWayFiles::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* datatype
 ) const

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

@@ -65,9 +65,9 @@ private:
 
   // private member functions
 
-    fileName getFilePath(word&, bool) const;
+    fileName getFilePath(const word&, bool) const;
 
-    void renameFilePath(fileName&,word&) const;
+    void renameFilePath(const fileName&, const word&) const;
 
 public:
 
@@ -92,24 +92,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const { FatalError<<"function not implemented !!! twoWayFiles::getData!!!\n" << abort(FatalError); }
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype = ""
         ) const;

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

@@ -138,8 +138,8 @@ double twoWayMPI::DEMVariableValue(word variablename)
 
 void twoWayMPI::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label /*step*/
 ) const
@@ -149,8 +149,8 @@ void twoWayMPI::getData
 
 void twoWayMPI::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     int ** const& field,
     label /*step*/
 ) const
@@ -160,8 +160,8 @@ void twoWayMPI::getData
 
 void twoWayMPI::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* datatype
 ) const

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

@@ -107,24 +107,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const;
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype
         ) const;

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

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

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

@@ -178,8 +178,8 @@ twoWayMany2Many::~twoWayMany2Many()
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 void twoWayMany2Many::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label /*step*/
 ) const
@@ -231,8 +231,8 @@ void twoWayMany2Many::getData
 
 void twoWayMany2Many::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     int ** const& field,
     label /*step*/
 ) const
@@ -242,8 +242,8 @@ void twoWayMany2Many::getData
 
 void twoWayMany2Many::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* /*datatype*/
 ) const

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

@@ -157,24 +157,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const;
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype
         ) const;

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

@@ -253,8 +253,8 @@ twoWayOne2One::~twoWayOne2One()
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 void twoWayOne2One::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     label /*step*/
 ) const
@@ -361,8 +361,8 @@ void twoWayOne2One::getData
 
 void twoWayOne2One::getData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     int ** const& field,
     label /*step*/
 ) const
@@ -383,8 +383,8 @@ void twoWayOne2One::getData
 
 void twoWayOne2One::giveData
 (
-    word name,
-    word type,
+    const word& name,
+    const word& type,
     double ** const& field,
     const char* datatype
 ) const
@@ -966,6 +966,10 @@ void twoWayOne2One::extractCollected(T*& src, T**& dst, int width) const
     }
 }
 
+int twoWayOne2One::getNumberOfParticles() const
+{
+    return particleCloud_.numberOfParticles();
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -152,24 +152,24 @@ public:
     // Member Functions
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             label step
         ) const;
 
         void getData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             int ** const& field,
             label step
         ) const;
 
         void giveData
         (
-            word name,
-            word type,
+            const word& name,
+            const word& type,
             double ** const& field,
             const char* datatype
         ) const;
@@ -208,6 +208,8 @@ public:
         template <typename T>
         void extractCollected(T*&, T**&, int width=1) const;
 
+        int getNumberOfParticles() const;
+
         scalar getCG() const { return lmp->force->cg(); }
 };
 

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

@@ -49,7 +49,8 @@ heatTransferGranConduction::heatTransferGranConduction
     totalHeatFlux_(0.0),
     QPartPartName_(propsDict_.lookupOrDefault<word>("QPartPartName","QPartPart")),
     QPartPart_
-    (   IOobject
+    (
+        IOobject
         (
             QPartPartName_,
             sm.mesh().time().timeName(),
@@ -72,7 +73,7 @@ heatTransferGranConduction::heatTransferGranConduction
         ),
         sm.mesh(),
         dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
-	"zeroGradient"
+        "zeroGradient"
     ),
     partThermCondField_
     (
@@ -86,7 +87,7 @@ heatTransferGranConduction::heatTransferGranConduction
         ),
         sm.mesh(),
         dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
-	"zeroGradient"
+        "zeroGradient"
     ),
     partTempField_(sm.mesh().lookupObject<volScalarField>("partTemp")),
     prescribedVoidfractionFieldName_(propsDict_.lookupOrDefault<word>("prescribedVoidfractionFieldName","voidfraction")),
@@ -94,11 +95,11 @@ heatTransferGranConduction::heatTransferGranConduction
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","conductiveHeatFlux")),
-    partHeatFlux_(NULL),
     typePartThermCond_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
-    partThermCond_(NULL)
+    partThermCondRegName_(typeName + "partThermCond")
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_,1);
+    particleCloud_.registerParticleProperty<double**>(partThermCondRegName_,1);
 
     if (typePartThermCond_[0] < 0.0)
     {
@@ -127,25 +128,13 @@ heatTransferGranConduction::heatTransferGranConduction
 
 heatTransferGranConduction::~heatTransferGranConduction()
 {
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
-    particleCloud_.dataExchangeM().destroy(partThermCond_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void heatTransferGranConduction::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
-    particleCloud_.dataExchangeM().allocateArray(partThermCond_,initVal,1);
-}
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void heatTransferGranConduction::calcEnergyContribution()
 {
-   // realloc the arrays
-    allocateMyArrays();
-
     calcPartEffThermCond();
 
     QPartPart_ = fvc::laplacian(partEffThermCondField_,partTempField_);
@@ -156,6 +145,7 @@ void heatTransferGranConduction::calcEnergyContribution()
     scalar voidfraction(1);
 
     totalHeatFlux_ = 0.0;
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
@@ -199,6 +189,8 @@ void heatTransferGranConduction::calcPartThermCond()
 {
     label cellI=0;
     label partType = 1;
+    double**& partThermCond_ = particleCloud_.getParticlePropertyRef<double**>(partThermCondRegName_);
+
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
@@ -227,7 +219,8 @@ void heatTransferGranConduction::calcPartThermCond()
 
 void heatTransferGranConduction::heatFlux(label index, scalar vol, scalar voidfraction, scalar QPartPart)
 {
-        partHeatFlux_[index][0] = vol * QPartPart / (1.0 - voidfraction) ;
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+    partHeatFlux_[index][0] = vol * QPartPart / (1.0 - voidfraction) ;
 }
 
 void heatTransferGranConduction::giveData()
@@ -238,6 +231,7 @@ void heatTransferGranConduction::giveData()
         Info << "total conductive particle-particle heat flux [W] (Eulerian) = " << totalHeatFlux_ << endl;
     }
 
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
     particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
 }
 

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

@@ -54,7 +54,7 @@ protected:
 
     scalar totalHeatFlux_;
 
-    word QPartPartName_;
+    const word QPartPartName_;
 
     volScalarField QPartPart_;
 
@@ -64,23 +64,19 @@ protected:
 
     const volScalarField& partTempField_;
 
-    word prescribedVoidfractionFieldName_;
+    const word prescribedVoidfractionFieldName_;
 
     const volScalarField& prescribedVoidfraction_;
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
-    word partHeatFluxName_;
-
-    mutable double **partHeatFlux_;
+    const word partHeatFluxName_;
 
     scalarList typePartThermCond_;
 
-    mutable double **partThermCond_;
-
-    void allocateMyArrays() const;
+    const word partThermCondRegName_;
 
     void calcPartEffThermCond();
 
@@ -118,9 +114,9 @@ public:
 
         void addEnergyCoefficient(volScalarField&) const {}
 
-	void calcEnergyContribution();
+        void calcEnergyContribution();
 
-	void postFlow();
+        void postFlow();
 
 };
 

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

@@ -143,17 +143,25 @@ heatTransferGunn::heatTransferGunn
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partTempName_(propsDict_.lookup("partTempName")),
-    partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","convectiveHeatFlux")),
-    partHeatFlux_(NULL),
-    partHeatFluxCoeff_(NULL),
-    partRe_(NULL),
-    partNu_(NULL),
+    partHeatFluxCoeffRegName_(typeName + "partHeatFluxCoeff"),
+    partReRegName_(typeName + "partRe"),
+    partNuRegName_(typeName + "partNu"),
     scaleDia_(1.),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     maxTypeCG_(typeCG_.size())
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_,1);
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_,1);
+    if (implicit_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partHeatFluxCoeffRegName_,1);
+    }
+    if(verbose_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partReRegName_,1);
+        particleCloud_.registerParticleProperty<double**>(partNuRegName_,1);
+    }
 
     if (propsDict_.found("NusseltScalingFactor"))
     {
@@ -226,41 +234,16 @@ heatTransferGunn::heatTransferGunn
 
 heatTransferGunn::~heatTransferGunn()
 {
-    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
-    particleCloud_.dataExchangeM().destroy(partRe_,1);
-    particleCloud_.dataExchangeM().destroy(partNu_,1);
-    if (implicit_)
-    {
-        particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1);
-    }
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void heatTransferGunn::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
-    particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
-    if(implicit_)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
-    }
-
-    if(verbose_)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
-    }
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void heatTransferGunn::calcEnergyContribution()
 {
-   // realloc the arrays
-    allocateMyArrays();
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     // reset Scalar field
     QPartFluid_.primitiveFieldRef() = 0.0;
@@ -387,6 +370,8 @@ void heatTransferGunn::calcEnergyContribution()
 
                 if(verbose_)
                 {
+                    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
+                    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuRegName_);
                     partRe_[index][0] = Rep;
                     partNu_[index][0] = Nup;
                 }
@@ -433,6 +418,7 @@ void heatTransferGunn::calcEnergyContribution()
 
     if(implicit_)
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
         QPartFluidCoeff_.primitiveFieldRef() = 0.0;
 
         particleCloud_.averagingM().setScalarSum
@@ -448,6 +434,8 @@ void heatTransferGunn::calcEnergyContribution()
 
     if(verbose_)
     {
+        double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
+        double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuRegName_);
         ReField_.primitiveFieldRef() = 0.0;
         NuField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
@@ -515,6 +503,8 @@ scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) con
 void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid, scalar cg3)
 {
     scalar hAs = h * As * cg3;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     if (particleCloud_.getParticleEffVolFactors())
     {
@@ -529,6 +519,7 @@ void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid,
     }
     else
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
         partHeatFluxCoeff_[index][0] = hAs;
     }
 }
@@ -541,6 +532,7 @@ void heatTransferGunn::giveData()
         reduce(totalHeatFlux_, sumOp<scalar>());
         Info << "total convective particle-fluid heat flux [W] = " << totalHeatFlux_ << endl;
     }
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
     particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
 }
 
@@ -552,6 +544,9 @@ void heatTransferGunn::postFlow()
         scalar Tfluid(0.0);
         scalar Tpart(0.0);
         interpolationCellPoint<scalar> TInterpolator_(tempField_);
+        double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+        double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
 
         totalHeatFlux_ = 0.0;
 
@@ -585,6 +580,7 @@ void heatTransferGunn::postFlow()
 
 void heatTransferGunn::partTempField()
 {
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
     partTempField_.primitiveFieldRef() = 0.0;
     particleCloud_.averagingM().resetWeightFields();
     particleCloud_.averagingM().setScalarAverage
@@ -607,6 +603,8 @@ void heatTransferGunn::initPartTemp()
 {
     label cellI = 0;
     scalar T = 0.0;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];

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

@@ -68,11 +68,11 @@ protected:
 
     scalar NusseltScalingFactor_;
 
-    word QPartFluidName_;
+    const word QPartFluidName_;
 
     volScalarField QPartFluid_;
 
-    word QPartFluidCoeffName_;
+    const word QPartFluidCoeffName_;
 
     volScalarField QPartFluidCoeff_;
 
@@ -90,37 +90,33 @@ protected:
 
     dimensionedScalar partTempAve_;
 
-    word tempFieldName_;
+    const word tempFieldName_;
 
     const volScalarField& tempField_;           // ref to temperature field
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;        // ref to voidfraction field
 
     scalar maxSource_;                          // max (limited) value of src field
 
-    word velFieldName_;
+    const word velFieldName_;
 
     const volVectorField& U_;
 
-    word densityFieldName_;
+    const word densityFieldName_;
 
     const volScalarField& rho_;
 
-    word partTempName_;
+    const word partTempName_;
 
-    mutable double **partTemp_;
+    const word partHeatFluxName_;
 
-    word partHeatFluxName_;
+    const word partHeatFluxCoeffRegName_;
 
-    mutable double **partHeatFlux_;
+    const word partReRegName_;
 
-    mutable double **partHeatFluxCoeff_;
-
-    mutable double **partRe_;
-
-    mutable double **partNu_;
+    const word partNuRegName_;
 
     mutable scalar scaleDia_;
 
@@ -128,8 +124,6 @@ protected:
 
     const label maxTypeCG_;
 
-    void allocateMyArrays() const;
-
     void partTempField();
 
     scalar Nusselt(scalar, scalar, scalar) const;

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

@@ -142,17 +142,25 @@ heatTransferRanzMarshall::heatTransferRanzMarshall
     densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partTempName_(propsDict_.lookup("partTempName")),
-    partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","convectiveHeatFlux")),
-    partHeatFlux_(NULL),
-    partHeatFluxCoeff_(NULL),
-    partRe_(NULL),
-    partNu_(NULL),
+    partHeatFluxCoeffRegName_(typeName + "partHeatFluxCoeff"),
+    partReRegName_(typeName + "partRe"),
+    partNuRegName_(typeName + "partNu"),
     scaleDia_(1.),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     maxTypeCG_(typeCG_.size())
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partTempName_,1);
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_,1);
+    if (implicit_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partHeatFluxCoeffRegName_,1);
+    }
+    if(verbose_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partReRegName_,1);
+        particleCloud_.registerParticleProperty<double**>(partNuRegName_,1);
+    }
 
     if (propsDict_.found("NusseltScalingFactor"))
     {
@@ -225,41 +233,16 @@ heatTransferRanzMarshall::heatTransferRanzMarshall
 
 heatTransferRanzMarshall::~heatTransferRanzMarshall()
 {
-    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
-    particleCloud_.dataExchangeM().destroy(partRe_,1);
-    particleCloud_.dataExchangeM().destroy(partNu_,1);
-    if (implicit_)
-    {
-        particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1);
-    }
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void heatTransferRanzMarshall::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
-    particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
-    if(implicit_)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
-    }
-
-    if(verbose_)
-    {
-        particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
-        particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
-    }
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void heatTransferRanzMarshall::calcEnergyContribution()
 {
-   // realloc the arrays
-    allocateMyArrays();
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     // reset Scalar field
     QPartFluid_.primitiveFieldRef() = 0.0;
@@ -387,6 +370,8 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
                 if(verbose_)
                 {
+                    double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
+                    double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuRegName_);
                     partRe_[index][0] = Rep;
                     partNu_[index][0] = Nup;
                 }
@@ -427,6 +412,7 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
     if(implicit_)
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
         QPartFluidCoeff_.primitiveFieldRef() = 0.0;
 
         particleCloud_.averagingM().setScalarSum
@@ -442,6 +428,8 @@ void heatTransferRanzMarshall::calcEnergyContribution()
 
     if(verbose_)
     {
+        double**& partRe_ = particleCloud_.getParticlePropertyRef<double**>(partReRegName_);
+        double**& partNu_ = particleCloud_.getParticlePropertyRef<double**>(partNuRegName_);
         ReField_.primitiveFieldRef() = 0.0;
         NuField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
@@ -503,8 +491,10 @@ scalar heatTransferRanzMarshall::Nusselt(scalar voidfraction, scalar Rep, scalar
 void heatTransferRanzMarshall::heatFlux(label index, scalar h, scalar As, scalar Tfluid, scalar cg3)
 {
     scalar hAs = h * As * cg3;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
-    if (particleCloud_.getParticleEffVolFactors()) 
+    if (particleCloud_.getParticleEffVolFactors())
     {
         scalar effVolFac = particleCloud_.particleEffVolFactor(index);
         hAs *= effVolFac;
@@ -517,6 +507,7 @@ void heatTransferRanzMarshall::heatFlux(label index, scalar h, scalar As, scalar
     }
     else
     {
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
         partHeatFluxCoeff_[index][0] = hAs;
     }
 }
@@ -529,6 +520,7 @@ void heatTransferRanzMarshall::giveData()
         reduce(totalHeatFlux_, sumOp<scalar>());
         Info << "total convective particle-fluid heat flux [W] = " << totalHeatFlux_ << endl;
     }
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
     particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
 }
 
@@ -540,6 +532,9 @@ void heatTransferRanzMarshall::postFlow()
         scalar Tfluid(0.0);
         scalar Tpart(0.0);
         interpolationCellPoint<scalar> TInterpolator_(tempField_);
+        double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+        double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+        double**& partHeatFluxCoeff_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxCoeffRegName_);
 
         totalHeatFlux_ = 0.0;
 
@@ -573,6 +568,7 @@ void heatTransferRanzMarshall::postFlow()
 
 void heatTransferRanzMarshall::partTempField()
 {
+        double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
         partTempField_.primitiveFieldRef() = 0.0;
         particleCloud_.averagingM().resetWeightFields();
         particleCloud_.averagingM().setScalarAverage
@@ -595,6 +591,8 @@ void heatTransferRanzMarshall::initPartTemp()
 {
     label cellI = 0;
     scalar T = 0.0;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];

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

@@ -67,11 +67,11 @@ protected:
 
     scalar NusseltScalingFactor_;
 
-    word QPartFluidName_;
+    const word QPartFluidName_;
 
     volScalarField QPartFluid_;
 
-    word QPartFluidCoeffName_;
+    const word QPartFluidCoeffName_;
 
     volScalarField QPartFluidCoeff_;
 
@@ -89,37 +89,33 @@ protected:
 
     dimensionedScalar partTempAve_;
 
-    word tempFieldName_;
+    const word tempFieldName_;
 
     const volScalarField& tempField_;            // ref to temperature field
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;    // ref to voidfraction field
 
     scalar maxSource_;                           // max (limited) value of src field
 
-    word velFieldName_;
+    const word velFieldName_;
 
     const volVectorField& U_;
 
-    word densityFieldName_;
+    const word densityFieldName_;
 
     const volScalarField& rho_;
 
-    word partTempName_;
+    const word partTempName_;
 
-    mutable double **partTemp_;
+    const word partHeatFluxName_;
 
-    word partHeatFluxName_;
+    const word partHeatFluxCoeffRegName_;
 
-    mutable double **partHeatFlux_;
+    const word partReRegName_;
 
-    mutable double **partHeatFluxCoeff_;
-
-    mutable double **partRe_;
-
-    mutable double **partNu_;
+    const word partNuRegName_;
 
     mutable scalar scaleDia_;
 
@@ -127,8 +123,6 @@ protected:
 
     const label maxTypeCG_;
 
-    void allocateMyArrays() const;
-
     void partTempField();
 
     scalar Nusselt(scalar, scalar, scalar) const;
@@ -163,11 +157,11 @@ public:
 
         void addEnergyContribution(volScalarField&) const;
 
-	void addEnergyCoefficient(volScalarField&) const;
+        void addEnergyCoefficient(volScalarField&) const;
 
-	void calcEnergyContribution();
+        void calcEnergyContribution();
 
-	void postFlow();
+        void postFlow();
 };
 
 

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

@@ -50,7 +50,6 @@ reactionHeat::reactionHeat
     mesh_(sm.mesh()),
     maxSource_(1e30),
     reactionHeatName_(propsDict_.lookupOrDefault<word>("reactionHeatName","reactionHeat")),
-    reactionHeat_(NULL),
     reactionHeatField_
     (
         IOobject
@@ -65,7 +64,7 @@ reactionHeat::reactionHeat
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0),0.0)
     )
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(reactionHeatName_,1);
 
     if(propsDict_.found("maxsource"))
     {
@@ -79,23 +78,15 @@ reactionHeat::reactionHeat
 
 reactionHeat::~reactionHeat()
 {
-    particleCloud_.dataExchangeM().destroy(reactionHeat_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void reactionHeat::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(reactionHeat_,initVal,1);
-}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 void reactionHeat::calcEnergyContribution()
 {
-   // realloc the arrays
-    allocateMyArrays();
+    double**& reactionHeat_ = particleCloud_.getParticlePropertyRef<double**>(reactionHeatName_);
 
     particleCloud_.dataExchangeM().getData(reactionHeatName_,"scalar-atom",reactionHeat_);
 

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

@@ -56,12 +56,8 @@ protected:
 
     word reactionHeatName_;
 
-    mutable double **reactionHeat_;
-
     volScalarField reactionHeatField_;
 
-    void allocateMyArrays() const;
-
 public:
 
     //- Runtime type information

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

@@ -54,7 +54,11 @@ FinesFields::FinesFields
     velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+#if OPENFOAM_VERSION_MAJOR < 5
+    voidfraction_(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_))),
+#else
     voidfraction_(sm.mesh().lookupObjectRef<volScalarField> (voidfractionFieldName_)),
+#endif
     UsFieldName_(propsDict_.lookupOrDefault<word>("granVelFieldName","Us")),
     UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
     pFieldName_(propsDict_.lookupOrDefault<word>("pFieldName","p")),
@@ -224,7 +228,7 @@ FinesFields::FinesFields
     g_("g",dimensionSet(0,1,-2,0,0),vector(0,0,-9.81)),
     alphaDynMax_(0.1),
     alphaMax_(propsDict_.lookupOrDefault<scalar>("alphaMax",0.95)),
-    alphaMinClog_(propsDict_.lookupOrDefault<scalar>("alphaMinClog",0.3)),
+    alphaMinClog_(propsDict_.lookupOrDefault<scalar>("alphaMinClog",0.1)),
     critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 0.05)),
     deltaT_(voidfraction_.mesh().time().deltaTValue()),
     depositionLength_(0.0),
@@ -412,7 +416,7 @@ void FinesFields::calcSource()
     {
         fKin = 0.0;
         fStick = 0.0;
-        if (clogKin_ && alphaP_[cellI] > alphaMinClog_) // no kinetic cloggig in dilute regions
+        if (clogKin_ && alphaP_[cellI] > alphaMinClog_) // no cloggig in dilute regions
         {
             // calculate everything in units auf dSauter
             critpore = nCrit_*dFine_.value()/dSauter_[cellI];
@@ -428,7 +432,7 @@ void FinesFields::calcSource()
             else if (fKin > 1.0) fKin = 1.0;
         }
 
-        if (clogStick_)
+        if (clogStick_ && alphaP_[cellI] > alphaMinClog_) // no cloggig in dilute regions
         {
             magU = mag(uReconstructed_()[cellI]); // use U reconstructed from phi to suppress oscillations at interfaces
            // fStick = 1.0 / ( 1.0 + magU/uBind_) * alphaP_[cellI] / 0.65;

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

@@ -74,8 +74,8 @@ KochHillRWDrag::KochHillRWDrag
     interpolation_(propsDict_.found("interpolation")),
     scale_(1.),
     randomTauE_(propsDict_.found("randomTauE")),
-    partTime_(NULL),
-    partUfluct_(NULL),
+    partTimeRegName_(typeName + "partTime"),
+    partUfluctRegName_(typeName + "partUfluct"),
     RanGen_(label(0))
 {
 
@@ -99,16 +99,8 @@ KochHillRWDrag::KochHillRWDrag
     if (propsDict_.found("rhoP"))
         rhoP_= readScalar(propsDict_.lookup("rhoP"));
 
-
-//    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-//    {
-        //allocate memory
-        particleCloud_.dataExchangeM().allocateArray(partTime_,0.,1);
-        particleCloud_.dataExchangeM().allocateArray(partUfluct_,0.,3);
-//    }
-
-    //Pout << "RW-TEST: maxNumberOfParticles() == " << particleCloud_.dataExchangeM().maxNumberOfParticles() << endl; // TEST-Output
-
+    particleCloud_.registerParticleProperty<double**>(partTimeRegName_,1,0.0,false);
+    particleCloud_.registerParticleProperty<double**>(partUfluctRegName_,3,0.0,false);
 }
 
 
@@ -116,18 +108,12 @@ KochHillRWDrag::KochHillRWDrag
 
 KochHillRWDrag::~KochHillRWDrag()
 {
-    particleCloud_.dataExchangeM().destroy(partTime_, 1);
-    particleCloud_.dataExchangeM().destroy(partUfluct_, 3);
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void KochHillRWDrag::setForce() const
 {
-
-    // realloc the arrays
-    reAllocArrays();
-
     if (scale_ > 1.0)
     {
         Info << "KochHillRW using scale = " << scale_ << endl;
@@ -189,6 +175,9 @@ void KochHillRWDrag::setForce() const
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
 
+    double**& partTime_ = particleCloud_.getParticlePropertyRef<double**>(partTimeRegName_);
+    double**& partUfluct_ = particleCloud_.getParticlePropertyRef<double**>(partUfluctRegName_);
+
     //Info << "RW-TEST: We are in setForce() at t = " << t << endl; // TEST-Output
 
     for (int index = 0; index<particleCloud_.numberOfParticles(); ++index)
@@ -382,16 +371,6 @@ void KochHillRWDrag::setForce() const
 }
 
 
-void KochHillRWDrag::reAllocArrays() const
-{
-    if (particleCloud_.numberOfParticlesChanged())
-    {
-        particleCloud_.dataExchangeM().allocateArray(partTime_,0.0,1);  // field/initVal/with/lenghtFromLigghts
-        particleCloud_.dataExchangeM().allocateArray(partUfluct_,0.0,3);
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -89,9 +89,9 @@ private:
 
     const bool randomTauE_;
 
-    mutable double **partTime_;                // Lagrangian array
+    const word partTimeRegName_;
 
-    mutable double **partUfluct_;              // Lagrangian array
+    const word partUfluctRegName_;
 
     mutable Random RanGen_;
 
@@ -117,8 +117,6 @@ public:
 
     // Member Functions
         void setForce() const;
-
-        void reAllocArrays() const;
 };
 
 

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

@@ -71,14 +71,10 @@ LaEuScalarTemp::LaEuScalarTemp
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     partTempName_(propsDict_.lookup("partTempName")),
-    partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookup("partHeatFluxName")),
-    partHeatFlux_(NULL),
     lambda_(readScalar(propsDict_.lookup("lambda"))),
     Cp_(readScalar(propsDict_.lookup("Cp")))
 {
-    allocateMyArrays();
-
     if (propsDict_.found("maxSource"))
     {
         maxSource_=readScalar(propsDict_.lookup ("maxSource"));
@@ -96,8 +92,10 @@ LaEuScalarTemp::LaEuScalarTemp
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
 
-
     particleCloud_.checkCG(false);
+
+    particleCloud_.registerParticleProperty<double**>(partTempName_);
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_);
 }
 
 
@@ -105,8 +103,6 @@ LaEuScalarTemp::LaEuScalarTemp
 
 LaEuScalarTemp::~LaEuScalarTemp()
 {
-    particleCloud_.dataExchangeM().destroy(partTemp_,1);
-    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -114,6 +110,9 @@ void LaEuScalarTemp::allocateMyArrays() const
 {
     // get memory for 2d arrays
     double initVal = 0.0;
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+
     particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
     particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
 }
@@ -128,6 +127,8 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
 {
     // realloc the arrays
     allocateMyArrays();
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
 
     // reset Scalar field
     EuField.primitiveFieldRef() = 0.0;

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

@@ -81,12 +81,8 @@ private:
 
     word partTempName_;
 
-    mutable double **partTemp_;                  // Lagrangian array
-
     word partHeatFluxName_;
 
-    mutable double **partHeatFlux_;              // Lagrangian array
-
     scalar lambda_;                              // fluid thermal conductivity [W/(m*K)]
 
     scalar Cp_;                                  // specific heat capacity [W*s/(kg*K)]

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

@@ -66,7 +66,9 @@ SchillerNaumannDrag::SchillerNaumannDrag
     propsDict_(dict.subDict(typeName + "Props")),
     verbose_(propsDict_.found("verbose")),
     velFieldName_(propsDict_.lookup("velFieldName")),
-    U_(sm.mesh().lookupObject<volVectorField> (velFieldName_))
+    U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
+    scaleDia_(1.),
+    scaleDrag_(1.)
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
@@ -85,7 +87,12 @@ SchillerNaumannDrag::SchillerNaumannDrag
     // read those switches defined above, if provided in dict
     forceSubM(0).readSwitches();
 
-    particleCloud_.checkCG(false);
+    particleCloud_.checkCG(true);
+
+    if (propsDict_.found("scale"))
+        scaleDia_ = scalar(readScalar(propsDict_.lookup("scale")));
+    if (propsDict_.found("scaleDrag"))
+        scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag")));
 }
 
 
@@ -99,6 +106,16 @@ SchillerNaumannDrag::~SchillerNaumannDrag()
 
 void SchillerNaumannDrag::setForce() const
 {
+    if (scaleDia_ > 1.0)
+    {
+        Info << "SchillerNaumann using scale = " << scaleDia_ << endl;
+    }
+    else if (particleCloud_.cg() > 1.0)
+    {
+        scaleDia_ = particleCloud_.cg();
+        Info << "SchillerNaumann using scale from liggghts cg = " << scaleDia_ << endl;
+    }
+
     #include "setupProbeModel.H"
 
     const volScalarField& nufField = forceSubM(0).nuField();
@@ -117,6 +134,7 @@ void SchillerNaumannDrag::setForce() const
                 vector Us = particleCloud_.velocity(index);
                 vector Ur = U_[cellI]-Us;
                 scalar ds = 2*particleCloud_.radius(index);
+                scalar ds_scaled = ds/scaleDia_;
                 scalar nuf = nufField[cellI];
                 scalar rho = rhoField[cellI];
                 scalar voidfraction = particleCloud_.voidfraction(index);
@@ -127,30 +145,33 @@ void SchillerNaumannDrag::setForce() const
                 if (magUr > 0)
                 {
                    // calc particle Re Nr
-                    Rep = ds*magUr/nuf;
+                    Rep = ds_scaled*magUr/nuf;
 
                     // calc fluid drag Coeff
                     Cd = max(0.44,24.0/Rep*(1.0+0.15*pow(Rep,0.687)));
 
                     // calc particle's drag
-                    drag = 0.125*Cd*rho*M_PI*ds*ds*magUr*Ur;
+                    drag = 0.125*Cd*rho*M_PI*ds*ds*scaleDia_*magUr*Ur*scaleDrag_;
 
                     if (modelType_=="B")
                         drag /= voidfraction;
                 }
 
-                if(verbose_ && index >100 && index <102)
+                if(verbose_ && index >=100 && index <102)
                 {
-                    Info << "index = " << index << endl;
-                    Info << "Us = " << Us << endl;
-                    Info << "Ur = " << Ur << endl;
-                    Info << "ds = " << ds << endl;
-                    Info << "rho = " << rho << endl;
-                    Info << "nuf = " << nuf << endl;
-                    Info << "voidfraction = " << voidfraction << endl;
-                    Info << "Rep = " << Rep << endl;
-                    Info << "Cd = " << Cd << endl;
-                    Info << "drag = " << drag << endl;
+                    Pout << "cellI = " << cellI << endl;
+                    Pout << "index = " << index << endl;
+                    Pout << "Ufluid = " << U_[cellI] << endl;
+                    Pout << "Us = " << Us << endl;
+                    Pout << "Ur = " << Ur << endl;
+                    Pout << "ds = " << ds << endl;
+                    Pout << "ds/scale = " << ds_scaled << endl;
+                    Pout << "rho = " << rho << endl;
+                    Pout << "nuf = " << nuf << endl;
+                    Pout << "voidfraction = " << voidfraction << endl;
+                    Pout << "Rep = " << Rep << endl;
+                    Pout << "Cd = " << Cd << endl;
+                    Pout << "drag = " << drag << endl;
                 }
 
                 //Set value fields and write the probe

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

@@ -66,6 +66,10 @@ private:
 
     const volVectorField& U_;
 
+    mutable scalar scaleDia_;
+
+    mutable scalar scaleDrag_;
+
 public:
 
     //- Runtime type information

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

@@ -54,8 +54,8 @@ dSauter::dSauter
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
     multiTypes_(false),
-    d2_(NULL),
-    d3_(NULL),
+    d2RegName_(typeName + "d2"),
+    d3RegName_(typeName + "d3"),
     maxTypeCG_(1),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
     d2Field_
@@ -101,9 +101,11 @@ dSauter::dSauter
         multiTypes_ = true;
         maxTypeCG_ = typeCG_.size();
     }
-    allocateMyArrays();
-    dSauter_.write();
 
+    particleCloud_.registerParticleProperty<double**>(d2RegName_,1);
+    particleCloud_.registerParticleProperty<double**>(d3RegName_,1);
+
+    dSauter_.write();
 
     // init force sub model
     setForceSubModels(propsDict_);
@@ -114,19 +116,10 @@ dSauter::dSauter
 
 dSauter::~dSauter()
 {
-    particleCloud_.dataExchangeM().destroy(d2_,1);
-    particleCloud_.dataExchangeM().destroy(d3_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
-void dSauter::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal = 0.0;
-    particleCloud_.dataExchangeM().allocateArray(d2_,initVal,1);  // field/initVal/with/lenghtFromLigghts
-    particleCloud_.dataExchangeM().allocateArray(d3_,initVal,1);
-}
 
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 
@@ -137,7 +130,8 @@ void dSauter::setForce() const
         Info << "dSauter using CG factor(s) = " << typeCG_ << endl;
     }
 
-    allocateMyArrays();
+    double**& d2_ = particleCloud_.getParticlePropertyRef<double**>(d2RegName_);
+    double**& d3_ = particleCloud_.getParticlePropertyRef<double**>(d3RegName_);
 
     label cellI = 0;
     label partType = 1;
@@ -151,11 +145,11 @@ void dSauter::setForce() const
         cellI = particleCloud_.cellIDs()[index][0];
         if (cellI >= 0)
         {
-            if (particleCloud_.getParticleEffVolFactors()) 
+            if (particleCloud_.getParticleEffVolFactors())
             {
                 effVolFac = particleCloud_.particleEffVolFactor(index);
             }
-            if (multiTypes_) 
+            if (multiTypes_)
             {
                 partType = particleCloud_.particleType(index);
                 if (partType > maxTypeCG_)

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

@@ -45,21 +45,19 @@ private:
 
     bool multiTypes_;
 
-    mutable double **d2_;
-    
-    mutable double **d3_;
+    const word d2RegName_;
+
+    const word d3RegName_;
 
     label maxTypeCG_;
 
     scalarList typeCG_;
-    
-    mutable volScalarField d2Field_;
-    
-    mutable volScalarField d3Field_;
-    
-    mutable volScalarField dSauter_;
 
-    void allocateMyArrays() const;
+    mutable volScalarField d2Field_;
+
+    mutable volScalarField d3Field_;
+
+    mutable volScalarField dSauter_;
 
 public:
 

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

@@ -53,7 +53,7 @@ granKineticEnergy::granKineticEnergy
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
-    vfluc_(NULL),
+    vflucRegName_(typeName + "vfluc_mag"),
     UsFieldName_(propsDict_.lookup("granVelFieldName")),
     UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
     granKineticEnergy_
@@ -70,7 +70,7 @@ granKineticEnergy::granKineticEnergy
         "zeroGradient"
     )
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(vflucRegName_,1);
     granKineticEnergy_.write();
 
 
@@ -82,21 +82,14 @@ granKineticEnergy::granKineticEnergy
 
 granKineticEnergy::~granKineticEnergy()
 {
-    particleCloud_.dataExchangeM().destroy(vfluc_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void granKineticEnergy::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal = 0.0;
-    particleCloud_.dataExchangeM().allocateArray(vfluc_,initVal,1);
-}
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 
 void granKineticEnergy::setForce() const
 {
-    allocateMyArrays();
+    double**& vfluc_ = particleCloud_.getParticlePropertyRef<double**>(vflucRegName_);
 
     label cellI = 0;
     vector velfluc(0,0,0);

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

@@ -43,16 +43,14 @@ private:
 
     dictionary propsDict_;
 
-    mutable double **vfluc_;
+    const word vflucRegName_;
 
-    word UsFieldName_;
+    const word UsFieldName_;
 
     const volVectorField& UsField_;
 
     mutable volScalarField granKineticEnergy_;
 
-    void allocateMyArrays() const;
-
 public:
 
     //- Runtime type information

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

@@ -94,8 +94,17 @@ particleCellVolume::particleCellVolume
     ),
     upperThreshold_(readScalar(propsDict_.lookup("upperThreshold"))),
     lowerThreshold_(readScalar(propsDict_.lookup("lowerThreshold"))),
-    verbose_(propsDict_.found("verbose"))
+    verbose_(propsDict_.found("verbose")),
+    writeToFile_(propsDict_.lookupOrDefault<Switch>("writeToFile",false)),
+    filePtr_()
 {
+    // create the path and output file
+    if(writeToFile_)
+    {
+        fileName path(particleCloud_.IOM().createTimeDir("postProcessing/particleCellVolume"));
+        filePtr_.set(new OFstream(path/"particleCellVolume.txt"));
+        filePtr_() << "# time | total particle volume in cells | total volume of cells with particles | average volume fraction | min(voidfraction) | max(voidfraction)" << endl;
+    }
 }
 
 
@@ -120,6 +129,8 @@ void particleCellVolume::setForce() const
 
         scalar fieldValue=-1;
         scalar cellVol=-1;
+        scalar minFieldVal=1e18;
+        scalar maxFieldVal=-1e18;
 
         forAll(field,cellI)
         {
@@ -129,6 +140,8 @@ void particleCellVolume::setForce() const
                 cellVol = mesh_.V()[cellI];
                 scalarField_[cellI] = (1-fieldValue) * cellVol;
                 scalarField2_[cellI] = cellVol;
+                minFieldVal = min(minFieldVal, fieldValue);
+                maxFieldVal = max(maxFieldVal, fieldValue);
             }
             else
             {
@@ -138,6 +151,8 @@ void particleCellVolume::setForce() const
         }
         scalarField_.ref() = gSum(scalarField_);
         scalarField2_.ref() = gSum(scalarField2_);
+        reduce(minFieldVal, minOp<scalar>());
+        reduce(maxFieldVal, maxOp<scalar>());
 
         if(verbose_)
         {
@@ -147,8 +162,20 @@ void particleCellVolume::setForce() const
                 << ", and > " << lowerThreshold_
                 << ",\n the total volume of cells holding particles = " << scalarField2_[0]
                 << ",\n this results in an average volume fraction of:" << scalarField_[0]/(scalarField2_[0]+SMALL)
+                << ",\n the min occurring " << scalarFieldName_ << " is:" << minFieldVal
+                << ",\n the max occurring " << scalarFieldName_ << " is:" << maxFieldVal
                 << endl;
         }
+
+        if(writeToFile_)
+        {
+            filePtr_()  << mesh_.time().value() << " "
+                        << scalarField_[0] << " "
+                        << scalarField2_[0] << " "
+                        << scalarField_[0]/(scalarField2_[0]+SMALL) << " "
+                        << minFieldVal << " "
+                        << maxFieldVal << endl;
+        }
     }// end if time >= startTime_
 }
 

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

@@ -74,6 +74,10 @@ private:
 
     const Switch verbose_;
 
+    const Switch writeToFile_;
+
+    mutable autoPtr<OFstream> filePtr_;
+
 public:
 
     //- Runtime type information

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

@@ -60,9 +60,9 @@ particleDeformation::particleDeformation
     partTypes_(propsDict_.lookupOrDefault<labelList>("partTypes",labelList(1,-1))),
     lowerBounds_(propsDict_.lookupOrDefault<scalarList>("lowerBounds",scalarList(1,-1.0))),
     upperBounds_(propsDict_.lookupOrDefault<scalarList>("upperBounds",scalarList(1,-1.0))),
-    partDeformations_(NULL)
+    partDeformationsName_("partDeformations")
 {
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(partDeformationsName_,1);
 
     // init force sub model
     setForceSubModels(propsDict_);
@@ -120,16 +120,9 @@ particleDeformation::particleDeformation
 
 particleDeformation::~particleDeformation()
 {
-    particleCloud_.dataExchangeM().destroy(partDeformations_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void particleDeformation::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal = 0.0;
-    particleCloud_.dataExchangeM().allocateArray(partDeformations_,initVal,1);
-}
 
 bool particleDeformation::defaultDeformCell(label cell) const
 {
@@ -145,8 +138,8 @@ void particleDeformation::setForce() const
         init();
         initialExec_ = false;
     }
-    // realloc the arrays
-    allocateMyArrays();
+
+    double**& partDeformations_ = particleCloud_.getParticlePropertyRef<double**>(partDeformationsName_);
 
     label cellI = 0;
     label partType = -1;
@@ -204,7 +197,7 @@ void particleDeformation::setForce() const
     }
 
     // give DEM data
-    particleCloud_.dataExchangeM().giveData("partDeformations","scalar-atom", partDeformations_);
+    particleCloud_.dataExchangeM().giveData(partDeformationsName_,"scalar-atom", partDeformations_);
 }
 
 void particleDeformation::init() const

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

@@ -55,12 +55,12 @@ private:
 
     mutable bool initialExec_;
 
-    word refFieldName_;
+    const word refFieldName_;
 
     mutable autoPtr<volScalarField> refField_;
 
     // default deformation in region
-    word defaultDeformCellsName_;
+    const word defaultDeformCellsName_;
 
     autoPtr<cellSet> defaultDeformCells_;
 
@@ -74,12 +74,10 @@ private:
 
     scalarList upperBounds_;
 
-    mutable double **partDeformations_;
+    const word partDeformationsName_;
 
     label getListIndex(label) const;
 
-    void allocateMyArrays() const;
-
     void init() const;
 
     bool defaultDeformCell(label) const;

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

@@ -18,7 +18,7 @@ SourceFiles
     pdCorrelation.C
 
 Contributing Author
-    2018    Paul Kieckhefen, TUHH 
+    2018    Paul Kieckhefen, TUHH
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
@@ -55,11 +55,11 @@ pdCorrelation::pdCorrelation
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
-    d_(nullptr),
-    p_(nullptr),
-    d2_(nullptr),
-    pd_(nullptr),
-    cg3_(nullptr),
+    dRegName_(typeName + "d"),
+    pRegName_(typeName + "p"),
+    d2RegName_(typeName + "d2"),
+    pdRegName_(typeName + "pd"),
+    cg3RegName_(typeName + "cg3"),
     dField_
     (   IOobject
         (
@@ -151,7 +151,11 @@ pdCorrelation::pdCorrelation
                   << abort(FatalError);
     }
 
-    allocateMyArrays();
+    particleCloud_.registerParticleProperty<double**>(dRegName_,1);
+    particleCloud_.registerParticleProperty<double**>(pRegName_,3);
+    particleCloud_.registerParticleProperty<double**>(d2RegName_,1);
+    particleCloud_.registerParticleProperty<double**>(pdRegName_,3);
+    particleCloud_.registerParticleProperty<double**>(cg3RegName_,1);
 
     dField_.write();
     pdField_.write();
@@ -165,24 +169,9 @@ pdCorrelation::pdCorrelation
 
 pdCorrelation::~pdCorrelation()
 {
-    particleCloud_.dataExchangeM().destroy(cg3_, 1);
-    particleCloud_.dataExchangeM().destroy(d_,  1);
-    particleCloud_.dataExchangeM().destroy(p_,  3);
-    particleCloud_.dataExchangeM().destroy(d2_, 1);
-    particleCloud_.dataExchangeM().destroy(pd_, 3);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void pdCorrelation::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal = 0.0;
-    particleCloud_.dataExchangeM().allocateArray(d_,  initVal, 1);
-    particleCloud_.dataExchangeM().allocateArray(p_,  initVal, 3);
-    particleCloud_.dataExchangeM().allocateArray(d2_, initVal, 1);
-    particleCloud_.dataExchangeM().allocateArray(pd_, initVal, 3);
-    particleCloud_.dataExchangeM().allocateArray(cg3_, initVal, 1);
-}
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 
 void pdCorrelation::setForce() const
@@ -191,7 +180,11 @@ void pdCorrelation::setForce() const
 
     if (runOnWriteOnly_ && !mesh.write()) return; // skip if it's not write time
 
-    allocateMyArrays();
+    double**& d_ = particleCloud_.getParticlePropertyRef<double**>(dRegName_);
+    double**& p_ = particleCloud_.getParticlePropertyRef<double**>(pRegName_);
+    double**& d2_ = particleCloud_.getParticlePropertyRef<double**>(d2RegName_);
+    double**& pd_ = particleCloud_.getParticlePropertyRef<double**>(pdRegName_);
+    double**& cg3_ = particleCloud_.getParticlePropertyRef<double**>(cg3RegName_);
 
     const Switch densityFromList
     (

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

@@ -45,11 +45,11 @@ private:
 
     dictionary propsDict_;
 
-    mutable double **d_;
-    mutable double **p_;
-    mutable double **d2_;
-    mutable double **pd_;
-    mutable double **cg3_;
+    const word dRegName_;
+    const word pRegName_;
+    const word d2RegName_;
+    const word pdRegName_;
+    const word cg3RegName_;
 
     mutable volScalarField dField_;
     mutable volVectorField pField_;
@@ -65,8 +65,6 @@ private:
     const Switch     CG_;
     const Switch     runOnWriteOnly_;
 
-    void allocateMyArrays() const;
-
 public:
 
     //- Runtime type information

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

@@ -1,6 +1,6 @@
 /*---------------------------------------------------------------------------*\
     CFDEMcoupling academic - Open Source CFD-DEM coupling
-    
+
     Contributing authors:
     Thomas Lichtenegger
     Copyright (C) 2015- Johannes Kepler University, Linz
@@ -83,10 +83,10 @@ potentialRelaxation::potentialRelaxation
     ignoreReg_(propsDict_.lookupOrDefault<bool>("ignoreRegion",false)),
     ignoreDirection_(propsDict_.lookupOrDefault<vector>("ignoreDirection",vector::zero)),
     ignorePoint_(propsDict_.lookupOrDefault<vector>("ignorePoint",vector::zero)),
-    vfluc_(NULL)
+    vflucName_("vfluc")
 {
-    allocateMyArrays();
-    
+    particleCloud_.registerParticleProperty<double**>(vflucName_,3);
+
     if(ignoreReg_)
     {
         if(mag(ignoreDirection_) < SMALL)
@@ -103,48 +103,39 @@ potentialRelaxation::potentialRelaxation
 
 potentialRelaxation::~potentialRelaxation()
 {
-    delete vfluc_;
 }
 
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void potentialRelaxation::allocateMyArrays() const
-{
-    // get memory for 2d arrays
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(vfluc_,initVal,3); 
-}
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void potentialRelaxation::setForce() const
 {
-
     relax(D0_,D1_);
-    
+
     volVectorField relaxStream = -fvc::grad(correctedField_);
-    
+
   // volVectorField relaxStream = DField_ * fvc::grad(voidfraction_ - voidfractionRec_);
-    
-     // realloc the arrays
-    allocateMyArrays();
-    
+
+    double**& vfluc_ = particleCloud_.getParticlePropertyRef<double**>(vflucName_);
+
     vector position(0,0,0);
     scalar voidfraction(0.0);
     vector flucU(0,0,0);
     label cellI=0;
-   
+
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> relaxStreamInterpolator_(relaxStream);
-    
+
     scalar dtDEM = particleCloud_.dataExchangeM().DEMts();
     scalar dtCFD = voidfraction_.mesh().time().deltaTValue();
-    
+
     // if DEM time step > CFD time step, scale velocity down
     scalar timeFac = 1.0;
     if (dtDEM > dtCFD) timeFac = dtCFD / dtDEM;
 
-    
+
     for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
@@ -164,18 +155,18 @@ void potentialRelaxation::setForce() const
                     {
                       position = particleCloud_.position(index);
                       voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
-                      flucU = relaxStreamInterpolator_.interpolate(position,cellI);      
+                      flucU = relaxStreamInterpolator_.interpolate(position,cellI);
                     }
                     else
                     {
                         voidfraction = voidfraction_[cellI];
                         flucU = relaxStream[cellI];
                     }
-                    
+
                     if (voidfraction > 1.0-SMALL) voidfraction = 1.0 - SMALL;
                     flucU /= (1-voidfraction);
-		    flucU *= timeFac;
-                    // write particle based data to global array 
+                    flucU *= timeFac;
+                    // write particle based data to global array
                     for(int i = 0; i < 3; i++)
                     {
                         vfluc_[index][i]=flucU[i];
@@ -183,13 +174,13 @@ void potentialRelaxation::setForce() const
                 }
             }
     }
-    
-    particleCloud_.dataExchangeM().giveData("vfluc","vector-atom", vfluc_);
-    
+
+    particleCloud_.dataExchangeM().giveData(vflucName_,"vector-atom", vfluc_);
+
     if (measureDiff_)
     {
         dimensionedScalar diff( fvc::domainIntegrate( sqr( voidfraction_ - voidfractionRec_ ) ) );
-        scalar t = particleCloud_.mesh().time().timeOutputValue(); 
+        scalar t = particleCloud_.mesh().time().timeOutputValue();
         recErrorFile_ << t << "\t" << diff.value() << endl;
     }
 }
@@ -198,12 +189,12 @@ void potentialRelaxation::relax(scalar D0, scalar D1) const
 {
     volScalarField src0 = voidfraction_ - voidfractionRec_;
     volScalarField src1 = voidfraction_ - voidfractionRec_;
-    
+
     forAll(src1, cellI)
     {
         if(src1[cellI] > 0.0) src1[cellI] = 0.0;
     }
-    
+
     solve
     (
         fvm::laplacian(correctedField_)

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

@@ -1,6 +1,6 @@
 /*---------------------------------------------------------------------------*\
     CFDEMcoupling academic - Open Source CFD-DEM coupling
-    
+
     Contributing authors:
     Thomas Lichtenegger
     Copyright (C) 2015- Johannes Kepler University, Linz
@@ -44,44 +44,42 @@ class potentialRelaxation
 {
 private:
     dictionary propsDict_;
-    
+
     bool interpolate_;
-    
+
     bool measureDiff_;
-    
+
     mutable OFstream recErrorFile_;
-    
-    word voidfractionFieldName_;
-    
+
+    const word voidfractionFieldName_;
+
     const volScalarField& voidfraction_;
-    
-    word voidfractionRecFieldName_;
-    
+
+    const word voidfractionRecFieldName_;
+
     const volScalarField& voidfractionRec_;
-    
+
     scalar critVoidfraction_;
-    
+
     scalar D0_;
-    
+
     scalar D1_;
-    
+
     mutable volScalarField correctedField_;
-    
+
     const scalar dt_;
-    
+
     // ignore particles in cells below plane given with ref point and normal vector
     // normal vector points towards region where fluctuations are permitted
-    
-    bool ignoreReg_;
-    
-    vector ignoreDirection_;
-    
-    vector ignorePoint_;
-    
-    mutable double **vfluc_;              // Lagrangian array
 
-    void allocateMyArrays() const;
-    
+    bool ignoreReg_;
+
+    vector ignoreDirection_;
+
+    vector ignorePoint_;
+
+    const word vflucName_;
+
     void relax(scalar, scalar) const;
 
 public:

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

@@ -0,0 +1,256 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "turbulentDispersion.H"
+#include "addToRunTimeSelectionTable.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(turbulentDispersion, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    turbulentDispersion,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+turbulentDispersion::turbulentDispersion
+(
+    const dictionary& dict,
+    cfdemCloud& sm,
+    word type
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(type + "Props")),
+    interpolate_(propsDict_.lookupOrDefault<bool>("interpolation", false)),
+    mesh_(sm.mesh()),
+    ignoreCellsName_(propsDict_.lookupOrDefault<word>("ignoreCellsName","none")),
+    ignoreCells_(),
+    existIgnoreCells_(true),
+    wallIndicatorField_
+    (   IOobject
+        (
+            "wallIndicator",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+    ),
+    minTurbKineticEnergy_(propsDict_.lookupOrDefault<scalar>("minTurbKineticEnergy", 0.0)),
+    turbKineticEnergyFieldName_(propsDict_.lookupOrDefault<word>("turbKineticEnergyFieldName","")),
+    turbKineticEnergy_(NULL),
+    existTurbKineticEnergyInObjReg_(false),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 0.9)),
+    ranGen_(clock::getTime()+pid())
+{
+    if (ignoreCellsName_ != "none")
+    {
+       ignoreCells_.set(new cellSet(particleCloud_.mesh(),ignoreCellsName_));
+       Info << type << ": ignoring fluctuations in cellSet " << ignoreCells_().name() <<
+        " with " << ignoreCells_().size() << " cells." << endl;
+    }
+    else existIgnoreCells_ = false;
+
+    // define a field to indicate if a cell is next to boundary
+     label cellI = -1;
+     forAll (mesh_.boundary(),patchI)
+     {
+         word patchName = mesh_.boundary()[patchI].name();
+         if (patchName.rfind("procB",0) == 0) continue;
+
+         forAll(mesh_.boundary()[patchI], faceI)
+         {
+             cellI = mesh_.boundary()[patchI].faceCells()[faceI];
+             wallIndicatorField_[cellI] = 1.0;
+         }
+     }
+
+    if (turbKineticEnergyFieldName_ != "")
+    {
+        existTurbKineticEnergyInObjReg_ = true;
+        volScalarField& k(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> (turbKineticEnergyFieldName_)));
+        turbKineticEnergy_ = &k;
+    }
+    else
+    {
+        turbKineticEnergy_ = new volScalarField
+        (
+            IOobject
+            (
+                "turbKinEnergy",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(0,2,-2,0,0), 0)
+        );
+    }
+
+    // make sure this is the last force model in list so that fluid velocity does not get overwritten
+    label numLastForceModel = sm.nrForceModels();
+    word lastForceModel = sm.forceModels()[numLastForceModel-1];
+    if (lastForceModel != "turbulentDispersion")
+    {
+        FatalError <<"Force model 'turbulentDispersion' needs to be last in list!\n" << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+turbulentDispersion::~turbulentDispersion()
+{
+    if (!existTurbKineticEnergyInObjReg_) delete turbKineticEnergy_;
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+bool turbulentDispersion::ignoreCell(label cell) const
+{
+    if (!existIgnoreCells_) return false;
+    else return ignoreCells_()[cell];
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void turbulentDispersion::setForce() const
+{
+    if (!existTurbKineticEnergyInObjReg_)
+    {
+        *turbKineticEnergy_ = particleCloud_.turbulence().k()();
+    }
+
+    label cellI = -1;
+    label patchID = -1;
+    label faceIGlobal = -1;
+    scalar flucProjection = 0.0;
+    scalar k = 0.0;
+    vector faceINormal = vector::zero;
+    vector flucU = vector::zero;
+    vector position = vector::zero;
+    word patchName("");
+
+    interpolationCellPoint<scalar> turbKineticEnergyInterpolator_(*turbKineticEnergy_);
+
+    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    {
+        cellI = particleCloud_.cellIDs()[index][0];
+        if (cellI > -1 && !ignoreCell(cellI))
+        {
+            // particles in dilute regions follow fluid without fluctuations
+
+            if (voidfraction_[cellI] < critVoidfraction_)
+            {
+                if (interpolate_)
+                {
+                    position = particleCloud_.position(index);
+                    k = turbKineticEnergyInterpolator_.interpolate(position,cellI);
+                }
+                else
+                {
+                    k = (*turbKineticEnergy_)[cellI];
+                }
+
+                if (k < minTurbKineticEnergy_) k = minTurbKineticEnergy_;
+
+                flucU=unitFlucDir()*Foam::sqrt(2.0*k);
+
+                // prevent particles being pushed through walls by regulating velocity fluctuations
+                // check if cell is adjacent to wall and remove corresponding components
+                if (wallIndicatorField_[cellI] > 0.5)
+                {
+                    const cell& faces = mesh_.cells()[cellI];
+                    forAll (faces, faceI)
+                    {
+                        faceIGlobal = faces[faceI];
+                        patchID = mesh_.boundaryMesh().whichPatch(faceIGlobal);
+                        if (patchID < 0) continue;
+                        patchName = mesh_.boundary()[patchID].name();
+
+                        if (patchName.rfind("procB",0) == 0) continue;
+
+                        faceINormal = mesh_.Sf()[faceIGlobal];
+                        faceINormal /= mag(faceINormal);
+                        flucProjection = faceINormal&flucU;
+                        if (flucProjection > 0.0) flucU -= flucProjection*faceINormal;
+                    }
+                }
+
+                for(int j=0;j<3;j++)
+                {
+                    particleCloud_.fluidVels()[index][j] += flucU[j];
+                }
+            }
+        }
+    }
+}
+
+vector turbulentDispersion::unitFlucDir() const
+{
+    // unit random vector
+    // algorithm according to:
+    // Marsaglia. "Choosing a point from the surface of a sphere." The Annals of Mathematical Statistics 43.2 (1972): 645-646.
+    scalar v1(0.0);
+    scalar v2(0.0);
+    scalar s(10.0);
+    scalar s2(0.0);
+    vector rvec(0,0,0);
+    while(s>1.0)
+    {
+        v1=2*(ranGen_.scalar01()-0.5);
+        v2=2*(ranGen_.scalar01()-0.5);
+        s=v1*v1+v2*v2;
+    }
+    s2=Foam::sqrt(1-s);
+    rvec[0]=2*v1*s2;
+    rvec[1]=2*v2*s2;
+    rvec[2]=1-2*s;
+    return rvec;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#ifndef turbulentDispersion_H
+#define turbulentDispersion_H
+
+#include "forceModel.H"
+#include "interpolationCellPoint.H"
+#include "Random.H"
+#include "autoPtr.H"
+#include "cellSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class turbulentDispersion Declaration
+\*---------------------------------------------------------------------------*/
+
+class turbulentDispersion
+:
+    public forceModel
+{
+protected:
+    dictionary propsDict_;
+
+    bool interpolate_;
+
+    const fvMesh& mesh_;
+
+    // ignore fluctuations in region
+    word ignoreCellsName_;
+
+    autoPtr<cellSet> ignoreCells_;
+
+    bool existIgnoreCells_;
+
+    mutable volScalarField wallIndicatorField_;
+
+    scalar minTurbKineticEnergy_;
+
+    word turbKineticEnergyFieldName_;
+
+    volScalarField* turbKineticEnergy_;
+
+    bool existTurbKineticEnergyInObjReg_;
+
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;
+
+    scalar critVoidfraction_;
+
+    virtual vector unitFlucDir() const;
+
+    mutable Random ranGen_;
+
+    bool ignoreCell(label) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("turbulentDispersion");
+
+
+    // Constructors
+
+        //- Construct from components
+        turbulentDispersion
+        (
+            const dictionary& dict,
+            cfdemCloud& sm,
+            word type = "turbulentDispersion"
+        );
+
+    // Destructor
+
+        ~turbulentDispersion();
+
+
+    // Member Functions
+        void setForce() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -69,16 +69,11 @@ virtualMassForce::virtualMassForce
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     phiFieldName_(propsDict_.lookup("phiFieldName")),
     phi_(sm.mesh().lookupObject<surfaceScalarField> (phiFieldName_)),
-    UrelOld_(NULL),
+    UrelOldRegName_(typeName + "UrelOld"),
     splitUrelCalculation_(propsDict_.lookupOrDefault<bool>("splitUrelCalculation",false)),
     Cadd_(0.5)
 {
-
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d array
-        particleCloud_.dataExchangeM().allocateArray(UrelOld_,NOTONCPU,3);
-    }
+    particleCloud_.registerParticleProperty<double**>(UrelOldRegName_,3,NOTONCPU,false);
 
     // init force sub model
     setForceSubModels(propsDict_);
@@ -117,15 +112,13 @@ virtualMassForce::virtualMassForce
 
 virtualMassForce::~virtualMassForce()
 {
-    particleCloud_.dataExchangeM().destroy(UrelOld_,3);
 }
 
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void virtualMassForce::setForce() const
 {
-    reAllocArrays();
+    double**& UrelOld_ = particleCloud_.getParticlePropertyRef<double**>(UrelOldRegName_);
 
     scalar dt = U_.mesh().time().deltaT().value();
 
@@ -233,17 +226,6 @@ void virtualMassForce::setForce() const
     }
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::virtualMassForce::reAllocArrays() const
-{
-    if(particleCloud_.numberOfParticlesChanged())
-    {
-        Pout << "virtualMassForce::reAllocArrays..." << endl;
-        particleCloud_.dataExchangeM().allocateArray(UrelOld_,NOTONCPU,3);
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -71,7 +71,7 @@ private:
 
     const surfaceScalarField& phi_;
 
-    mutable double **UrelOld_;
+    const word UrelOldRegName_;
 
     const bool splitUrelCalculation_;      //indicator to split calculation of Urel between CFDEM and LIGGGHTS
                                            //requires the integration fix to take dv/dt into account!
@@ -100,8 +100,6 @@ public:
 
     // Member Functions
         void setForce() const;
-
-        void reAllocArrays() const;
 };
 
 

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

@@ -69,7 +69,9 @@ volWeightedAverage::volWeightedAverage
     vectorFieldNames_(propsDict_.lookup("vectorFieldNames")),
     upperThreshold_(readScalar(propsDict_.lookup("upperThreshold"))),
     lowerThreshold_(readScalar(propsDict_.lookup("lowerThreshold"))),
-    verbose_(propsDict_.found("verbose"))
+    verbose_(propsDict_.found("verbose")),
+    writeToFile_(propsDict_.lookupOrDefault<Switch>("writeToFile", false)),
+    filePtr_()
 {
     // create vol weighted average scalar fields
     scalarFields_.setSize(scalarFieldNames_.size());
@@ -124,6 +126,13 @@ volWeightedAverage::volWeightedAverage
             )
         );
     }
+
+    // create the path and output file
+    if(writeToFile_)
+    {
+        fileName path(particleCloud_.IOM().createTimeDir("postProcessing/volWeightedAverage"));
+        filePtr_.set(new OFstream(path/"volWeightedAverage.txt"));
+    }
 }
 
 
@@ -141,6 +150,8 @@ void volWeightedAverage::setForce() const
     {
         if(verbose_) Info << "volWeightedAverage.C - setForce()" << endl;
 
+        if(writeToFile_) filePtr_() << mesh_.time().value() << " ";
+
         for (int i=0;i < scalarFieldNames_.size(); i++)
         {
             // get reference to actual field
@@ -171,7 +182,7 @@ void volWeightedAverage::setForce() const
             MPI_Allreduce(&totVol, &totVol_all, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
             integralValue = gSum(scalarFields_[i]);
             volWeightedAverage = integralValue / (totVol_all+SMALL);
-            scalarFields_[i].ref() = volWeightedAverage;
+            scalarFields_[i].primitiveFieldRef() = volWeightedAverage;
 
             if(verbose_)
             {
@@ -184,6 +195,8 @@ void volWeightedAverage::setForce() const
                      << ",\n considering cells where the field < " << upperThreshold_
                      << ", and > " << lowerThreshold_ << endl;
             }
+
+            if(writeToFile_) filePtr_() << volWeightedAverage << " ";
         }
 
         for (int i=0;i < vectorFieldNames_.size(); i++)
@@ -225,7 +238,11 @@ void volWeightedAverage::setForce() const
                      << ",\n considering cells where the mag(field) < " << upperThreshold_
                      << ", and > " << lowerThreshold_ << endl;
             }
+
+            if(writeToFile_) filePtr_() << volWeightedAverage << " ";
         }
+
+        if(writeToFile_) filePtr_() << endl;
     }// end if time >= startTime_
 }
 

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

@@ -76,6 +76,10 @@ private:
 
     const Switch verbose_;
 
+    const Switch writeToFile_;
+
+    mutable autoPtr<OFstream> filePtr_;
+
 public:
 
     //- Runtime type information

src/lagrangian/cfdemParticle/subModels/regionModel/allRegion/allRegion.C

@@ -78,7 +78,6 @@ allRegion::~allRegion()
 
 void allRegion::defineRegion() const
 {
-    reAllocArrays();
     // do nothing
 }
 

src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C

@@ -46,15 +46,6 @@ defineRunTimeSelectionTable(regionModel, dictionary);
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void regionModel::reAllocArrays() const
-{
-    if(particleCloud_.numberOfParticlesChanged())
-    {
-        // get arrays of new length
-        particleCloud_.dataExchangeM().allocateArray(inRegion_,1.,1);
-        particleCloud_.dataExchangeM().allocateArray(outRegion_,1.,1);
-    }
-}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -66,16 +57,10 @@ regionModel::regionModel
 )
 :
     dict_(dict),
-    particleCloud_(sm),
-    inRegion_(NULL),
-    outRegion_(NULL)
+    particleCloud_(sm)
 {
-    if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
-    {
-        // get memory for 2d arrays
-        particleCloud_.dataExchangeM().allocateArray(inRegion_,1.,1);
-        particleCloud_.dataExchangeM().allocateArray(outRegion_,1.,1);
-    }
+    particleCloud_.registerParticleProperty<double**>("inRegion",1,1.0,false);
+    particleCloud_.registerParticleProperty<double**>("outRegion",1,1.0,false);
 }
 
 
@@ -83,8 +68,6 @@ regionModel::regionModel
 
 regionModel::~regionModel()
 {
-    particleCloud_.dataExchangeM().destroy(inRegion_,1);
-    particleCloud_.dataExchangeM().destroy(outRegion_,1);
 }
 
 

src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.H

@@ -60,10 +60,6 @@ protected:
 
         cfdemCloud& particleCloud_;
 
-        mutable double **inRegion_;
-
-        mutable double **outRegion_;
-
 public:
 
     friend class voidFractionModel;
@@ -117,13 +113,6 @@ public:
 
         virtual void resetVolFields(volVectorField&) const = 0;
 
-        void reAllocArrays() const;
-
-    // Access
-        inline double ** const& inRegion()const{ return inRegion_; };
-
-        inline double ** const& outRegion()const { return outRegion_; };
-
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -55,14 +55,17 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     typeKs_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
-    partKs_(NULL)
+    partKsRegName_(typeName + "partKs")
 {
     if (typeKs_[0] < 0.0)
     {
         FatalError << "ZehnerSchluenderThermCond: provide list of thermal conductivities." << abort(FatalError);
     }
 
-    if (typeKs_.size() > 1) allocateMyArrays();
+    if (typeKs_.size() > 1)
+    {
+        particleCloud_.registerParticleProperty<double**>(partKsRegName_,1);
+    }
     else
     {
         partKsField_ *= typeKs_[0];
@@ -74,14 +77,6 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
 
 ZehnerSchluenderThermCond::~ZehnerSchluenderThermCond()
 {
-    if (typeKs_.size() > 1) particleCloud_.dataExchangeM().destroy(partKs_,1);
-}
-
-
-void ZehnerSchluenderThermCond::allocateMyArrays() const
-{
-    double initVal=0.0;
-    particleCloud_.dataExchangeM().allocateArray(partKs_,initVal,1);
 }
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -133,7 +128,7 @@ void ZehnerSchluenderThermCond::calcPartKsField() const
         FatalError << "ZehnerSchluenderThermCond needs data for more than one type, but types are not communicated." << abort(FatalError);
     }
 
-    allocateMyArrays();
+    double**& partKs_ = particleCloud_.getParticlePropertyRef<double**>(partKsRegName_);
     label cellI=0;
     label partType = 0;
     for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)

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

@@ -54,19 +54,17 @@ private:
 
     dictionary propsDict_;
 
-    word partKsFieldName_;
+    const word partKsFieldName_;
 
     volScalarField& partKsField_;
 
-    word voidfractionFieldName_;
+    const word voidfractionFieldName_;
 
     const volScalarField& voidfraction_;
 
     scalarList typeKs_;
 
-    mutable double **partKs_;
-
-    void allocateMyArrays() const;
+    const word partKsRegName_;
 
     void calcPartKsField() const;
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C

@@ -89,8 +89,6 @@ GaussVoidFraction::~GaussVoidFraction()
 
 void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV)
 {
-    reAllocArrays();
-
     voidfractionNext_.ref()=1;
 
     scalar radius(-1);
@@ -103,12 +101,12 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
         //if(mask[index][0])
         //{
             //reset
-            for(int subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+            for(int subcell=0;subcell<cellsPerParticle()[index][0];subcell++)
             {
                 particleWeights[index][subcell]=0;
                 particleVolumes[index][subcell]=0;
             }
-            cellsPerParticle_[index][0]=1;
+            cellsPerParticle()[index][0]=1;
             particleV[index][0]=0;
 
             //collecting data
@@ -142,7 +140,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
                 }
                 else if (hashSetLength > 0)
                 {
-                    cellsPerParticle_[index][0]=hashSetLength;
+                    cellsPerParticle()[index][0]=hashSetLength;
 
                     //making sure that the cell containing the center is the first subcell
                     particleCloud_.cellIDs()[index][0]=particleCenterCellID;
@@ -200,7 +198,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
     //bringing eulerian field to particle array
     for(label index=0; index< particleCloud_.numberOfParticles(); index++)
     {
-        for(label subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+        for(label subcell=0;subcell<cellsPerParticle()[index][0];subcell++)
         {
             label cellID = particleCloud_.cellIDs()[index][subcell];
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C

@@ -91,8 +91,6 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 {
     const boundBox& globalBb = particleCloud_.mesh().bounds();
 
-    reAllocArrays();
-
     voidfractionNext_.ref() = 1.0;
 
     for (int index=0; index < particleCloud_.numberOfParticles(); index++)
@@ -100,13 +98,13 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
         //if(mask[index][0])
         //{
             //reset
-            for (int subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
+            for (int subcell=0; subcell < cellsPerParticle()[index][0]; subcell++)
             {
                 particleWeights[index][subcell] = 0.0;
                 particleVolumes[index][subcell] = 0.0;
             }
 
-            cellsPerParticle_[index][0]=1;
+            cellsPerParticle()[index][0]=1;
             particleV[index][0]=0;
 
             //collecting data
@@ -266,7 +264,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
                 }
                 else if (hashSetLength > 0)
                 {
-                    cellsPerParticle_[index][0]=hashSetLength;
+                    cellsPerParticle()[index][0]=hashSetLength;
                     hashSett.erase(particleCenterCellID);
 
                     for (label i=0; i < hashSetLength-1; i++)
@@ -281,7 +279,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
     for (label index=0; index < particleCloud_.numberOfParticles(); index++)
     {
-        for (label subcell=0; subcell < cellsPerParticle_[index][0]; subcell++)
+        for (label subcell=0; subcell < cellsPerParticle()[index][0]; subcell++)
         {
             label cellID = particleCloud_.cellIDs()[index][subcell];
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C

@@ -88,8 +88,6 @@ bigParticleVoidFraction::~bigParticleVoidFraction()
 
 void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV)
 {
-    reAllocArrays();
-
     voidfractionNext_.ref()=1;
 
     scalar radius(-1);
@@ -102,12 +100,12 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
         //if(mask[index][0])
         //{
             //reset
-            for(int subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+            for(int subcell=0;subcell<cellsPerParticle()[index][0];subcell++)
             {
                 particleWeights[index][subcell]=0;
                 particleVolumes[index][subcell]=0;
             }
-            cellsPerParticle_[index][0]=1;
+            cellsPerParticle()[index][0]=1;
             particleV[index][0]=0;
 
             //collecting data
@@ -137,7 +135,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
                 }
                 else if (hashSetLength > 0)
                 {
-                    cellsPerParticle_[index][0]=hashSetLength;
+                    cellsPerParticle()[index][0]=hashSetLength;
 
                     //making sure that the cell containing the center is the first subcell
                     particleCloud_.cellIDs()[index][0]=particleCenterCellID;
@@ -190,7 +188,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
     //bringing eulerian field to particle array
     for(label index=0; index< particleCloud_.numberOfParticles(); index++)
     {
-        for(label subcell=0;subcell<cellsPerParticle_[index][0];subcell++)
+        for(label subcell=0;subcell<cellsPerParticle()[index][0];subcell++)
         {
             label cellID = particleCloud_.cellIDs()[index][subcell];
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C

@@ -80,8 +80,6 @@ centreVoidFraction::~centreVoidFraction()
 
 void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV)
 {
-    reAllocArrays();
-
     scalar radius(-1);
     scalar volume(0);
     scalar cellVol(0);
@@ -93,7 +91,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
         //{
             // reset
             particleWeights[index][0]=0;
-            cellsPerParticle_[index][0]=1;
+            cellsPerParticle()[index][0]=1;
 
             label cellI = particleCloud_.cellIDs()[index][0];
 
@@ -128,9 +126,9 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
                 if (index==0)
                 {
                     Info << "centre cellI = " << cellI << endl;
-                    Info << "cellsPerParticle_=" << cellsPerParticle_[index][0] << endl;
+                    Info << "cellsPerParticle =" << cellsPerParticle()[index][0] << endl;
 
-                    for(int i=0;i<cellsPerParticle_[index][0];i++)
+                    for(int i=0;i<cellsPerParticle()[index][0];i++)
                     {
                        if(i==0)Info << "cellids, voidfractions, particleWeights, : \n";
                        Info << particleCloud_.cellIDs()[index][i] << " ," << endl;