merge PFM/master > TUe/master

.circleci/config.yml

@@ -21,18 +21,18 @@ jobs:
           
       - run:
           name: Make project and user dir
-          command: mkdir -p /root/CFDEM/CFDEMcoupling && mkdir -p /root/CFDEM/-4.1
+          command: mkdir -p /root/CFDEM/CFDEMcoupling && mkdir -p /root/CFDEM/-6
 
       - checkout:
           path: /root/CFDEM/CFDEMcoupling
 
       - run:
           name: Add OpenFOAM package repository
-          command: sudo apt-get install -y software-properties-common wget apt-transport-https && sudo add-apt-repository http://dl.openfoam.org/ubuntu && sudo sh -c "wget -O - http://dl.openfoam.org/gpg.key | apt-key add -"
+          command: sudo apt-get install -y software-properties-common wget apt-transport-https && sudo sh -c "wget -O - http://dl.openfoam.org/gpg.key | apt-key add -" && sudo add-apt-repository http://dl.openfoam.org/ubuntu
 
       - run:
-          name: Install OpenFOAM 4.1
+          name: Install OpenFOAM 6
-          command: sudo apt-get update && sudo apt-get -y install openfoam4
+          command: sudo apt-get update && sudo apt-get -y install openfoam6
 
       - run:
           name: Clone LIGGGHTS repository
@@ -42,7 +42,7 @@ jobs:
           name: Build LIGGGHTS
           command: >
             shopt -s expand_aliases &&
-            source /opt/openfoam4/etc/bashrc &&
+            source /opt/openfoam6/etc/bashrc &&
             source /root/CFDEM/CFDEMcoupling/etc/bashrc &&
             bash /root/CFDEM/CFDEMcoupling/etc/compileLIGGGHTS.sh
           no_output_timeout: 30m
@@ -51,7 +51,7 @@ jobs:
           name: Build CFDEMcoupling library
           command: >
             shopt -s expand_aliases &&
-            source /opt/openfoam4/etc/bashrc &&
+            source /opt/openfoam6/etc/bashrc &&
             source /root/CFDEM/CFDEMcoupling/etc/bashrc &&
             bash /root/CFDEM/CFDEMcoupling/etc/compileCFDEMcoupling_src.sh
           
@@ -59,7 +59,7 @@ jobs:
           name: Build CFDEMcoupling solvers
           command: >
             shopt -s expand_aliases &&
-            source /opt/openfoam4/etc/bashrc &&
+            source /opt/openfoam6/etc/bashrc &&
             source /root/CFDEM/CFDEMcoupling/etc/bashrc &&
             bash /root/CFDEM/CFDEMcoupling/etc/compileCFDEMcoupling_sol.sh
 
@@ -67,6 +67,6 @@ jobs:
           name: Build CFDEMcoupling utilities
           command: >
             shopt -s expand_aliases &&
-            source /opt/openfoam4/etc/bashrc &&
+            source /opt/openfoam6/etc/bashrc &&
             source /root/CFDEM/CFDEMcoupling/etc/bashrc &&
             bash /root/CFDEM/CFDEMcoupling/etc/compileCFDEMcoupling_uti.sh

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

@@ -1,6 +1,10 @@
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
+
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
+     $(PFLAGS) \
     -I$(CFDEM_OFVERSION_DIR) \
     -ImultiphaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels \
@@ -13,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/cfdemSolverMultiphase.C

@@ -46,6 +46,11 @@ Description
 
 int main(int argc, char *argv[])
 {
+
+#if OPENFOAM_VERSION_MAJOR >= 6
+    FatalError << "cfdemSolverMultiphase requires OpenFOAM 4.x or 5.x to work properly" << exit(FatalError);
+#endif
+
     #include "postProcess.H"
     #include "setRootCase.H"
     #include "createTime.H"

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/Make/options

@@ -1,10 +1,15 @@
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
+
 EXE_INC = \
+     $(PFLAGS) \
     -IalphaContactAngle \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -Wno-deprecated-copy
 
 LIB_LIBS = \
     -linterfaceProperties \

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C

@@ -686,8 +686,13 @@ void Foam::multiphaseMixture::solveAlphas
             alphaPhiCorr,
             zeroField(),
             zeroField(),
+#if OPENFOAM_VERSION_MAJOR < 6
             1,
             0,
+#else
+            oneField(),
+            zeroField(),
+#endif
             true
         );
 

applications/solvers/cfdemSolverMultiphaseScalar/EEqn.H

@@ -1,7 +1,6 @@
         // get mixture properties
         Cp  = mixture.Cp();
         kf  = mixture.kf(); 
-        keff=particleCloud.thermCondM().thermCond(); 
 
         // get scalar source from DEM        
         particleCloud.energyContributions(Qsource);
@@ -11,7 +10,7 @@
         (
             rho*Cp*(fvm::ddt(voidfraction,T)
             + fvm::div(phi,T))
-	 - fvm::laplacian(keff*voidfraction,T)
+            - fvm::laplacian(thCond*voidfraction,T)
             ==
             Qsource + fvm::Sp(QCoeff,T)
         );

applications/solvers/cfdemSolverMultiphaseScalar/createFields.H

@@ -176,18 +176,19 @@ volScalarField kf
     mixture.kf()
 );
 
-volScalarField keff
+volScalarField thCond
 (
     IOobject
     (
-        "keff",
+        "thCond",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+    dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+    "zeroGradient"
 );
 
 Info<< "Reading/creating concentration fields\n" << endl;

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/EEqn.H

@@ -6,9 +6,6 @@
     particleCloud.energyContributions(Qsource);
     particleCloud.energyCoefficients(QCoeff);
 
-    //thDiff=particleCloud.thermCondM().thermDiff();
-    thCond=particleCloud.thermCondM().thermCond();
-
     addSource = fvc::ddt(rhoeps, K) + fvc::div(phi, K)
       + (
             he.name() == "e"
@@ -35,7 +32,6 @@
       - Qsource
       - fvm::Sp(QCoeff/Cpv, he)
      // thermal conduction of the fluid with effective conductivity
-     // - fvm::laplacian(rhoeps*thDiff,he)
       - fvm::laplacian(voidfraction*thCond/Cpv,he)
       // + particle-fluid energy transfer due to work
       // + fluid energy dissipation due to shearing
@@ -54,9 +50,6 @@
 
     thermo.correct();
 
-    Info << "Qsource" << max(Qsource).value() << " " << min(Qsource).value() << endl;
-    Info << "QCoeff" << max(QCoeff).value() << " " << min(QCoeff).value() << endl;
-    Info << "Cpv" << max(Cpv).value() << " " << min(Cpv).value() << endl;
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
 
     particleCloud.clockM().start(31,"energySolve");

applications/solvers/cfdemSolverRhoPimple/Make/options

@@ -1,5 +1,7 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
 PFLAGS+= -Dcompre
 
 EXE_INC = \
@@ -15,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/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C

@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
         particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
-        if(hasEvolved)
+        if(hasEvolved && smoothenForces)
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
         }
@@ -113,7 +113,11 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().start(26,"Flow");
 
+#if OPENFOAM_VERSION_MAJOR < 6
         if (pimple.nCorrPIMPLE() <= 1)
+#else
+        if (pimple.nCorrPimple() <= 1)
+#endif
         {
             #include "rhoEqn.H"
         }

applications/solvers/cfdemSolverRhoPimple/createFields.H

@@ -58,10 +58,11 @@ Info<< "Reading thermophysical properties\n" << endl;
            "addSource",
            runTime.timeName(),
            mesh,
-           IOobject::MUST_READ,
+           IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     );
 
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
@@ -94,21 +95,6 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
     );
 
-  /*  Info<< "\nCreating thermal diffusivity field\n" << endl;
-    volScalarField thDiff
-    (
-        IOobject
-        (
-            "thDiff",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
-    );
-  */
     Info<< "\nCreating thermal conductivity field\n" << endl;
     volScalarField thCond
     (
@@ -117,11 +103,12 @@ Info<< "Reading thermophysical properties\n" << endl;
             "thCond",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
     );
 
     Info<< "\nCreating heat capacity field\n" << endl;
@@ -190,6 +177,17 @@ Info<< "Reading thermophysical properties\n" << endl;
         )
     );
 
+    bool smoothenForces
+    (
+        pimple.dict().lookupOrDefault<bool>
+        (
+            "smoothenForces",
+            false
+        )
+    );
+    if (smoothenForces) Info << "Smoothening implicit particle forces.\n" << endl;
+    else Info << "Not smoothening implicit particle forces.\n" << endl;
+
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (

applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H

@@ -6,7 +6,6 @@ volScalarField& he = thermo.he();
 particleCloud.energyContributions(Qsource);
 particleCloud.energyCoefficients(QCoeff);
 
-thCond=particleCloud.thermCondM().thermCond();
 Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
 // correct source for the thermodynamic reference temperature

applications/solvers/cfdemSolverRhoPimpleChem/Make/options

@@ -1,5 +1,7 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
 PFLAGS+= -Dcompre
 
 EXE_INC = \
@@ -19,19 +21,14 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam \
-
+    -Wno-deprecated-copy
-
 
 
 EXE_LIBS = \
@@ -48,9 +45,6 @@ EXE_LIBS = \
     -l$(CFDEM_LIB_COMP_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS) \
-    -lliquidProperties \
-    -lliquidMixtureProperties \
-    -lthermophysicalFunctions \
     -lreactionThermophysicalModels \
     -lchemistryModel \
     -lradiationModels \

applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H

@@ -13,7 +13,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 {
     combustion->correct();
+#if OPENFOAM_VERSION_MAJOR < 5
     dQ = combustion->dQ();
+#else
+    Qdot = combustion->Qdot();
+#endif
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
@@ -72,4 +76,5 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         Y[inertIndex].max(0.0);
     }
 }
+
 particleCloud.clockM().stop("Y");

applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C

@@ -30,7 +30,12 @@ Description
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
+#if OPENFOAM_VERSION_MAJOR < 6
 #include "rhoCombustionModel.H"
+#else
+#include "rhoReactionThermo.H"
+#include "CombustionModel.H"
+#endif
 #include "bound.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
@@ -115,7 +120,11 @@ int main(int argc, char *argv[])
 
        particleCloud.clockM().start(26,"Flow");
 
+#if OPENFOAM_VERSION_MAJOR < 6
         if (pimple.nCorrPIMPLE() <= 1)
+#else
+        if (pimple.nCorrPimple() <= 1)
+#endif
         {
             #include "rhoEqn.H"
         }

applications/solvers/cfdemSolverRhoPimpleChem/createFields.H

@@ -1,20 +1,28 @@
     //  thermodynamics, chemistry
 
+#if OPENFOAM_VERSION_MAJOR < 6
     Info<< "Creating combustion model\n" << endl;
-
     autoPtr<combustionModels::rhoCombustionModel> combustion
     (
         combustionModels::rhoCombustionModel::New(mesh)
     );
-
     rhoReactionThermo& thermo = combustion->thermo();
+#else
+    Info<< "Reading thermophysical properties\n" << endl;
+    autoPtr<rhoReactionThermo> pThermo(rhoReactionThermo::New(mesh));
+    rhoReactionThermo& thermo = pThermo();
+#endif
     thermo.validate(args.executable(), "h", "e");
 
     basicSpecieMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
 
     // read molecular weight
+#if OPENFOAM_VERSION_MAJOR < 6
     volScalarField W(composition.W());
+#else
+    volScalarField W(thermo.W());
+#endif
 
     bool propagateInertSpecie = true;
 
@@ -127,11 +135,12 @@
             "thCond",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
     );
 
     Info<< "\nCreating heat capacity field\n" << endl;
@@ -197,6 +206,14 @@
         )
     );
 
+#if OPENFOAM_VERSION_MAJOR >= 6
+    Info<< "Creating combustion model\n" << endl;
+    autoPtr<CombustionModel<rhoReactionThermo>> combustion
+    (
+        CombustionModel<rhoReactionThermo>::New(thermo, turbulence())
+    );
+#endif
+
     Info<< "Creating field dpdt\n" << endl;
     volScalarField dpdt
     (
@@ -213,6 +230,7 @@
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
 
+#if OPENFOAM_VERSION_MAJOR < 5
     volScalarField dQ
     (
         IOobject
@@ -226,6 +244,21 @@
         mesh,
         dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
+#else
+    volScalarField Qdot
+    (
+        IOobject
+        (
+            "Qdot",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
+    );
+#endif
 
     Info<< "\nReading momentum exchange field Ksl\n" << endl;
     volScalarField Ksl

applications/solvers/cfdemSolverRhoSimple/Make/files

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

applications/solvers/rStatAnalysis/Make/files

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

applications/solvers/rStatAnalysis/Make/options

@@ -0,0 +1,27 @@
+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 \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/rStatAnalysis/rStatAnalysis.C

@@ -0,0 +1,62 @@
+/*---------------------------------------------------------------------------*\
+    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
+    rStatAnalysis
+
+Description
+    Creates and analyzes a recurrence statistics
+
+\*---------------------------------------------------------------------------*/
+
+#include "recBase.H"
+#include "recStatAnalysis.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+
+
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nAnalyzing recurrence statistics\n" << endl;
+
+    recurrenceBase.recStatA().init();
+    recurrenceBase.recStatA().statistics();
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverBase/Make/files

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

applications/solvers/rcfdemSolverBase/Make/options

@@ -0,0 +1,28 @@
+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 \
+    -Wno-deprecated-copy
+
+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/rcfdemSolverBase/createFields.H

@@ -0,0 +1,125 @@
+    // dummy fields
+    Info << "\nCreating dummy pressure and density fields\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0)
+    );
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+    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
+    );
+
+    // 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 phi\n" << endl;
+    surfaceScalarField phiRec
+    (
+        IOobject
+        (
+            "phiRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(URec*voidfractionRec) & mesh.Sf()
+    );
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );

applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    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
+    cfdemSolverRecurrence
+
+Description
+    Solves a transport equation for a passive scalar on a two-phase solution
+    Test-bed for a solver based on recurrence statistics
+
+Rules
+    Solution data to compute the recurrence statistics from, needs to
+        reside in $CASE_ROOT/dataBase
+    Time step data in dataBase needs to be evenly spaced in time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+
+#include "cfdemCloud.H"
+#include "clockModel.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<cfdemCloud> 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();
+
+    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");
+
+
+        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        {
+            recurrenceBase.updateRecFields();
+            #include "readFields.H"
+            recTimeIndex++;
+        }
+
+        particleCloud.clockM().start(27,"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/rcfdemSolverBase/readFields.H

@@ -0,0 +1,4 @@
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolVectorField("U",URec);
+recurrenceBase.recM().exportVolVectorField("Us",UsRec);
+recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);

applications/solvers/rcfdemSolverCoupledHeattransfer/Make/files

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

applications/solvers/rcfdemSolverCoupledHeattransfer/Make/options

@@ -0,0 +1,28 @@
+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 \
+    -Wno-deprecated-copy
+
+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/rcfdemSolverCoupledHeattransfer/TEqImp.H

@@ -0,0 +1,35 @@
+    volScalarField rhoeps = rhoRec*voidfractionRec;
+
+    particleCloud.energyContributions(Qsource);
+
+    particleCloud.energyCoefficients(QCoeff);
+
+    //K = 0.5*magSqr(URec);
+
+    addSource = fvc::div(phiRec/fvc::interpolate(rhoRec), pRec);
+    // main contribution due to gas expansion, not due to transport of kinetic energy
+    // fvc::ddt(rhoeps, K) + fvc::div(phiRec, K)
+
+    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
+    );
+
+    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/rcfdemSolverCoupledHeattransfer/createFields.H

@@ -0,0 +1,230 @@
+    // dummy fields
+    Info << "\nCreating dummy pressure field\n" << endl;
+    volScalarField pRec
+    (
+        IOobject
+        (
+            "pRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-2,0,0,0,0), 0.0)
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+
+    volScalarField rhoRec
+    (
+        IOobject
+        (
+            "rhoRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    volVectorField URec
+    (
+        IOobject
+        (
+            "URec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0, 1, -1, 0, 0), vector::zero)
+    );
+
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+    );
+
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0, 1, -1, 0, 0), vector::zero)
+    );
+
+    // 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
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,0,-1,0,0,0,0), 0.0)
+    );
+    phiRec.write();
+
+
+    Info << "Creating dummy turbulence model\n" << endl;
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );
+
+
+    const IOdictionary& transportProps = mesh.lookupObject<IOdictionary>("transportProperties");
+    dimensionedScalar molMass(transportProps.lookup("molM"));
+    // need to scale R down with 1e3 because return value of RR in g, not kg
+    dimensionedScalar R("R",dimensionSet(0,2,-2,-1,0,0,0),Foam::constant::thermodynamic::RR / (1e3*molMass.value()));
+    Info << "specific gas constant R = " << R << endl;
+
+    dimensionedScalar Cv(transportProps.lookup("Cv"));
+
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+           "addSource",
+           runTime.timeName(),
+           mesh,
+           IOobject::READ_IF_PRESENT,
+           IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+
+   // place to put weight functions
+    IOdictionary weightDict
+    (
+        IOobject
+        (
+            "weightDict",
+            runTime.constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    );
+    weightDict.add("weights",scalarList(1,1.0));

applications/solvers/rcfdemSolverCoupledHeattransfer/rcfdemSolverCoupledHeattransfer.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+    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 "fvOptions.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+#include "recPath.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);
+    #include "updateFields.H"
+    #include "updateRho.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    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();
+    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 "updateRho.H"
+        #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/rcfdemSolverCoupledHeattransfer/updateFields.H

@@ -0,0 +1,38 @@
+// get current weights for various databases
+
+// A: triggered over current value of boundary field
+// word boundaryName = "inlet";
+// label myinlet = mesh.boundary().findPatchID(boundaryName);
+// label startIndex = mesh.boundary()[boundaryName].start();
+
+// B: explicitly define weights
+
+scalarList wList(weightDict.lookupOrDefault("weights",scalarList(1,0.0)));
+
+recurrenceBase.recP().updateIntervalWeights(wList);
+
+// is it neccessary to extend recurrence path?
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolScalarField("p",pRec);
+recurrenceBase.recM().exportVolVectorField("Us",UsRec);
+recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);
+
+Info << "current database weights: = " << wList << endl;
+Info << "current database: " << recurrenceBase.recM().currDataBase() << endl;
+for(int i=0;i<wList.size();i++)
+{
+    scalar w = wList[i];
+    if (recurrenceBase.recM().currDataBase() == i) w -= 1.0;
+    phiRec += w*recurrenceBase.recM().exportSurfaceScalarFieldAve("phi",i)();
+}
+
+{
+    volScalarField& NuField(const_cast<volScalarField&>(mesh.lookupObject<volScalarField> ("NuField")));
+    recurrenceBase.recM().exportVolScalarField("NuField",NuField);
+}
+

applications/solvers/rcfdemSolverCoupledHeattransfer/updateRho.H

@@ -0,0 +1 @@
+rhoRec = pRec / (T * R);

applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H

@@ -6,9 +6,6 @@
     particleCloud.energyContributions(Qsource);
     particleCloud.energyCoefficients(QCoeff);
 
-    //thDiff=particleCloud.thermCondM().thermDiff();
-    thCond=particleCloud.thermCondM().thermCond();
-
     addSource =
         (
             he.name() == "e"
@@ -16,7 +13,7 @@
             fvc::div(phi, K) +
             fvc::div
             (
-                fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U),
+                fvc::absolute(phi/fvc::interpolate(rho), voidfractionRec*U),
                 p,
                 "div(phiv,p)"
             )
@@ -25,9 +22,6 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
-    // correct source for the thermodynamic reference temperature
-    dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
-    Qsource += QCoeff*Tref;
 
     fvScalarMatrix EEqn
     (
@@ -35,12 +29,12 @@
       + addSource
       - Qsource
       - fvm::Sp(QCoeff/Cpv, he)
-      - fvm::laplacian(voidfraction*thCond/Cpv,he)
+ //     - fvm::laplacian(voidfractionRec*kf/Cpv,he)
+      - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
 
-
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/rcfdemSolverRhoSteadyPimple/Make/files

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

applications/solvers/rcfdemSolverRhoSteadyPimple/Make/options

@@ -1,5 +1,7 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
 PFLAGS+= -Dcompre
 
 EXE_INC = \
@@ -15,9 +17,12 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -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)\
+    -lrecurrence \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \

applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H

@@ -4,7 +4,7 @@ particleCloud.otherForces(fOther);
 tmp<fvVectorMatrix> tUEqn
 (
     fvm::div(phi, U)
-  + particleCloud.divVoidfractionTau(U, voidfraction)
+  + particleCloud.divVoidfractionTau(U, voidfractionRec)
   + fvm::Sp(Ksl,U)
   - fOther
  ==
@@ -18,13 +18,16 @@ fvOptions.constrain(UEqn);
 
 if (modelType=="B" || modelType=="Bfull")
 {
-    solve(UEqn == -fvc::grad(p)+ Ksl*Us);
+    solve(UEqn == -fvc::grad(p)+ Ksl*UsRec);
 }
 else
 {
-    solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us);
+    solve(UEqn == -voidfractionRec*fvc::grad(p)+ Ksl*UsRec);
 }
 
+//U.relax();
+#include "limitU.H"
+
 fvOptions.correct(U);
 
 K = 0.5*magSqr(U);

applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H

@@ -51,6 +51,19 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh
     );
 
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        voidfraction
+    );
+
     volScalarField addSource
     (
         IOobject
@@ -58,10 +71,11 @@ Info<< "Reading thermophysical properties\n" << endl;
            "addSource",
            runTime.timeName(),
            mesh,
-           IOobject::MUST_READ,
+           IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     );
 
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
@@ -102,11 +116,12 @@ Info<< "Reading thermophysical properties\n" << endl;
             "thCond",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
     );
 
     Info<< "\nCreating heat capacity field\n" << endl;
@@ -158,7 +173,7 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar::lookupOrDefault
         (
             "rhoMax",
-            simple.dict(),
+            pimple.dict(),
             dimDensity,
             GREAT
         )
@@ -169,12 +184,45 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar::lookupOrDefault
         (
             "rhoMin",
-            simple.dict(),
+            pimple.dict(),
             dimDensity,
             0
         )
     );
 
+    dimensionedScalar pMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMax",
+            pimple.dict(),
+            dimPressure,
+            GREAT
+        )
+    );
+
+    dimensionedScalar pMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMin",
+            pimple.dict(),
+            dimPressure,
+            -GREAT
+        )
+    );
+
+    dimensionedScalar UMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "UMax",
+            pimple.dict(),
+            dimVelocity,
+            -1.0
+        )
+    );
+
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
@@ -189,7 +237,7 @@ Info<< "Reading thermophysical properties\n" << endl;
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
-    setRefCell(p, simple.dict(), pRefCell, pRefValue);
+    setRefCell(p, pimple.dict(), pRefCell, pRefValue);
 
     mesh.setFluxRequired(p.name());
 
@@ -217,11 +265,11 @@ Info<< "Reading thermophysical properties\n" << endl;
             "Ksl",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
-        //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
+        dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 0.0)
     );
 
 
@@ -239,4 +287,20 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh
     );
 
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        Us
+    );
+
+    
+    dimensionedScalar kf("0", dimensionSet(1, 1, -3, -1, 0, 0, 0), 0.026);
+
 //===============================

applications/solvers/rcfdemSolverRhoSteadyPimple/limitP.H

@@ -0,0 +1,2 @@
+p = max(p, pMin);
+p = min(p, pMax);

applications/solvers/rcfdemSolverRhoSteadyPimple/limitU.H

@@ -0,0 +1,11 @@
+if (UMax.value() > 0)
+{
+    forAll(U,cellI)
+    {
+        scalar mU(mag(U[cellI]));
+        if (mU > UMax.value())
+        {
+            U[cellI] *= UMax.value() / mU;
+        }
+    }
+}

applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H

@@ -7,14 +7,15 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
 if (modelType=="A")
 {
-    rhorAUf *= fvc::interpolate(voidfraction);
+    rhorAUf *= fvc::interpolate(voidfractionRec);
 }
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+//tUEqn.clear();
 
-surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
+surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*UsRec)& mesh.Sf());
 
 
-if (simple.transonic())
+if (pimple.transonic())
 {
 //  transonic version not implemented yet
 }
@@ -34,7 +35,7 @@ else
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, rhoeps, U, phi, rhorAUf);
 
-    while (simple.correctNonOrthogonal())
+    while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -49,7 +50,7 @@ else
 
         pEqn.solve();
 
-        if (simple.finalNonOrthogonalIter())
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
@@ -59,6 +60,8 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
+#include "limitP.H"
+
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
 rho = max(rho, rhoMin);
@@ -69,13 +72,15 @@ Info<< "rho max/min : " << max(rho).value()
 
 if (modelType=="A")
 {
-    U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);
+    U = HbyA - rAU*(voidfractionRec*fvc::grad(p)-Ksl*UsRec);
 }
 else
 {
-    U = HbyA - rAU*(fvc::grad(p)-Ksl*Us);
+    U = HbyA - rAU*(fvc::grad(p)-Ksl*UsRec);
 }
 
+#include "limitU.H"
+
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);

applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C

@@ -17,23 +17,31 @@ License
     Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
 
 Application
-    cfdemSolverRhoSimple
+    rcfdemSolverRhoSteadyPimple
 
 Description
-    Steady-state solver for turbulent flow of compressible fluids based on
+    Transient (DEM) + steady-state (CFD) solver for compressible flow using the
-    rhoSimpleFoam where functionality for CFD-DEM coupling has been added.
+    flexible PIMPLE (PISO-SIMPLE) algorithm. Particle-motion is obtained from
-
+    a recurrence process.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
+    where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
-#include "simpleControl.H"
+#include "pimpleControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+#include "recPath.H"
+
 #include "cfdemCloudEnergy.H"
 #include "implicitCouple.H"
 #include "clockModel.H"
@@ -61,18 +69,37 @@ int main(int argc, char *argv[])
     #include "createFvOptions.H"
 
     // create cfdemCloud
-    #include "readGravitationalAcceleration.H"
+    //#include "readGravitationalAcceleration.H"
-    cfdemCloudEnergy particleCloud(mesh);
+    cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
     #include "checkModelType.H"
+    recBase recurrenceBase(mesh);
+    #include "updateFields.H"
 
     turbulence->validate();
+    //#include "compressibleCourantNo.H"
+    //#include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+    label recTimeIndex = 0;
+    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
+    scalar startTime = runTime.startTime().value();
+
+    const IOdictionary& couplingProps = particleCloud.couplingProperties();
+    label nEveryFlow(couplingProps.lookupOrDefault<label>("nEveryFlow",1));
+    Info << "Solving flow equations every " << nEveryFlow << " steps.\n" << endl;
+    label stepcounter = 0;
+
     Info<< "\nStarting time loop\n" << endl;
 
-    while (simple.loop())
+    while (runTime.run())
     {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
         particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -81,6 +108,9 @@ int main(int argc, char *argv[])
         particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
+        //voidfraction = voidfractionRec;
+        //Us = UsRec;
+
         if(hasEvolved)
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
@@ -101,25 +131,56 @@ int main(int argc, char *argv[])
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
+        volScalarField rhoeps("rhoeps",rho*voidfractionRec);
+        if (stepcounter%nEveryFlow==0)
+        {
+            while (pimple.loop())
+            {
+                // if needed, perform drag update here
+#if OPENFOAM_VERSION_MAJOR < 6
+                if (pimple.nCorrPIMPLE() <= 1)
+#else
+                if (pimple.nCorrPimple() <= 1)
+#endif
+                {
+                    #include "rhoEqn.H"
+                }
 
-        volScalarField rhoeps("rhoeps",rho*voidfraction);
+                // --- Pressure-velocity PIMPLE corrector loop
-        // Pressure-velocity SIMPLE corrector
 
                 #include "UEqn.H"
-
-
-        // besides this pEqn, OF offers a "simple consistent"-option
-        #include "pEqn.H"
-        rhoeps=rho*voidfraction;
-
                 #include "EEqn.H"
 
-        turbulence->correct();
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    // besides this pEqn, OF offers a "pimple consistent"-option
+                    #include "pEqn.H"
+                    rhoeps=rho*voidfractionRec;
+                }
 
-        particleCloud.clockM().start(32,"postFlow");
+                if (pimple.turbCorr())
-        if(hasEvolved) particleCloud.postFlow();
+                {
+                    turbulence->correct();
+                }
+            }
+        }
+        stepcounter++;
+        particleCloud.clockM().stop("Flow");
+
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
         particleCloud.clockM().stop("postFlow");
 
+        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");
+
         runTime.write();
 
 
@@ -127,7 +188,7 @@ int main(int argc, char *argv[])
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
 
-        particleCloud.clockM().stop("Flow");
+
         particleCloud.clockM().stop("Global");
     }
 

applications/solvers/rcfdemSolverRhoSteadyPimple/updateFields.H

@@ -0,0 +1,8 @@
+// is it neccessary to extend recurrence path?
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolVectorField("Us",UsRec);

applications/solvers/recSolverTurbTransport/Make/files

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

applications/solvers/recSolverTurbTransport/Make/options

@@ -0,0 +1,27 @@
+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/recSolverTurbTransport/TEq.H

@@ -0,0 +1,17 @@
+
+volScalarField alphaEff("alphaEff", turbulence->nu()/Sc + dU2/Sct);
+
+TEqn =
+(
+    fvm::ddt(T)
+    + fvm::div(phiRec, T)
+    - fvm::laplacian(alphaEff, T)
+    ==
+    fvOptions(T)
+);
+
+TEqn.relax(relaxCoeff);
+
+fvOptions.constrain(TEqn);
+
+TEqn.solve();

applications/solvers/recSolverTurbTransport/createFields.H

@@ -0,0 +1,174 @@
+    // dummy fields
+    Info<< "\nCreating dummy pressure and density fields\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0)
+    );
+    
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+    
+    // recurrence fields
+    Info<< "\nCreating recurrence fields.\n" << endl;
+
+    volVectorField URec
+    (
+        IOobject
+        (
+            "URec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField U2Rec
+    (
+        IOobject
+        (
+            "U2Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // calculated fields
+    Info<< "\nCreating fields subject to calculation\n" << endl;
+
+    volScalarField delta
+    (
+        IOobject
+        (
+            "delta",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+       mesh,
+       dimensionedScalar("delta", dimLength, 0.0)
+    );
+
+    delta.primitiveFieldRef()=pow(mesh.V(),1.0/3.0);
+    delta.write();
+
+
+    Info<< "\ncreating dU2\n" << endl;
+
+    volScalarField dU2
+    (
+        IOobject
+        (
+            "dU2",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sqrt(0.5*mag(U2Rec - magSqr(URec)))*delta*0.094
+     );
+
+    forAll(dU2, cellI)
+    {
+        if (U2Rec[cellI]-magSqr(URec[cellI]) < 0.0)
+        {
+            dU2[cellI] = 0.0;
+        }
+    }
+
+    dU2.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) & mesh.Sf()
+ );
+
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );
+  
+    dimensionedScalar Sc("Sc", dimless, laminarTransport);
+    dimensionedScalar Sct("Sct", dimless, laminarTransport);
+  
+ // create concentration field
+ Info<< "Creating scalar transport field\n" << endl;
+
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    
+    fvScalarMatrix TEqn(T, dimless*dimVolume/(dimTime));
+
+    scalar relaxCoeff(0.0);
+
+ Info<< "reading clockProperties\n" << endl;
+
+    IOdictionary clockProperties
+    (
+            IOobject
+            (
+                "clockProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+    );
+
+    autoPtr<clockModel> myClock
+    (
+            clockModel::New
+            (
+                clockProperties,
+                mesh.time()
+            )
+    );
+

applications/solvers/recSolverTurbTransport/readFields.H

@@ -0,0 +1,14 @@
+
+recurrenceBase.recM().exportVolScalarField("U2Mean",U2Rec);
+recurrenceBase.recM().exportVolVectorField("UMean",URec);
+phiRec=linearInterpolate(URec) & mesh.Sf();
+
+dU2=sqrt(0.5*mag(U2Rec - magSqr(URec)))*delta*0.094;
+
+forAll(dU2, cellI)
+{
+    if (U2Rec[cellI]-magSqr(URec[cellI]) < 0.0)
+    {
+        dU2[cellI] = 0.0;
+    }
+}

applications/solvers/recSolverTurbTransport/recSolverTurbTransport.C

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    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
+    Turbulent Transport Solver Recurrence
+
+Description
+    Solves a transport equation for a passive scalar on a single-phase solution
+    for a solver based on recurrence statistics
+
+Rules
+	Solution data to compute the recurrence statistics from, needs to 
+		reside in $CASE_ROOT/dataBase
+	Time step data in dataBase needs to be evenly spaced in time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "clockModel.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"
+  
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+    
+    label recTimeIndex(0);
+    scalar recTimeStep_=recurrenceBase.recM().recTimeStep();
+    
+    while (runTime.run())
+    {
+
+        myClock().start(1,"Global");
+
+        runTime++;
+        
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        myClock().start(2,"fieldUpdate");
+        
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+	    Info << "Updating fields at run time " << runTime.timeOutputValue()
+	        << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            recurrenceBase.updateRecFields();
+            #include "readFields.H"
+            recTimeIndex++;
+        }
+
+        myClock().stop("fieldUpdate");
+
+        myClock().start(3,"speciesEqn");
+        #include "TEq.H"
+        myClock().stop("speciesEqn");
+
+        runTime.write();
+        
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+        
+        myClock().stop("Global");
+
+    }
+    
+    myClock().evalPar();
+    myClock().normHist();
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rtfmSolverSpecies/Make/files

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

applications/solvers/rtfmSolverSpecies/Make/options

@@ -0,0 +1,28 @@
+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 \
+    -Wno-deprecated-copy
+
+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/rtfmSolverSpecies/TEq.H

@@ -0,0 +1,14 @@
+    TEqn =
+    (
+        fvm::ddt(alpha2Rec, T)
+        + fvm::div(phi2Rec, T)
+        - fvm::laplacian(alpha2Rec*turbulence->nu(), T)
+        ==
+        fvOptions(alpha2Rec, T) // no fvOptions support yet
+    );
+
+    TEqn.relax(relaxCoeff);
+
+    fvOptions.constrain(TEqn); // no fvOptions support yet
+
+    TEqn.solve();

applications/solvers/rtfmSolverSpecies/createFields.H

@@ -0,0 +1,135 @@
+    // dummy fields
+    Info << "\nCreating dummy pressure and density fields\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0)
+    );
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+    volVectorField U1Rec
+    (
+        IOobject
+        (
+            "U1Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField alpha1Rec
+    (
+        IOobject
+        (
+            "alpha1Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volVectorField U2Rec
+    (
+        IOobject
+        (
+            "U2Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // calculated fields
+    Info << "\nCreating fields subject to calculation\n" << endl;
+    volScalarField alpha2Rec
+    (
+        IOobject
+        (
+            "alpha2Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        1-alpha1Rec
+    );
+
+
+    // write fields for t=t_start
+    alpha2Rec.write();
+//===============================
+
+
+    Info << "Calculating face flux field phi\n" << endl;
+    surfaceScalarField phi2Rec
+    (
+        IOobject
+        (
+            "phi2Rec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(U2Rec*alpha2Rec) & mesh.Sf()
+    );
+    phi2Rec.write();
+
+    singlePhaseTransportModel laminarTransport(U2Rec, phi2Rec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U2Rec, phi2Rec, laminarTransport)
+    );
+
+    // transport stuff
+
+    // create concentration field
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    fvScalarMatrix TEqn(T, dimless*dimVolume/(dimTime));
+
+    T.write();
+
+    scalar relaxCoeff(0.0);

applications/solvers/rtfmSolverSpecies/readFields.H

@@ -0,0 +1,6 @@
+recurrenceBase.recM().exportVolScalarField("alpha.air",alpha1Rec);
+alpha2Rec=1-alpha1Rec;
+recurrenceBase.recM().exportVolVectorField("U.air",U1Rec);
+recurrenceBase.recM().exportVolVectorField("U.water",U2Rec);
+recurrenceBase.recM().exportSurfaceScalarField("phi.water",phi2Rec);
+

applications/solvers/rtfmSolverSpecies/rtfmSolverSpecies.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    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
+    cfdemSolverRecurrence
+
+Description
+    Solves a transport equation for a passive scalar on a two-phase solution
+    Test-bed for a solver based on recurrence statistics
+
+Rules
+    Solution data to compute the recurrence statistics from, needs to
+        reside in $CASE_ROOT/dataBase
+    Time step data in dataBase needs to be evenly spaced in time
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+
+#include "cfdemCloud.H"
+#include "clockModel.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<cfdemCloud> particleCloud(mesh);
+    recBase recurrenceBase(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+
+    label recTimeIndex(0);
+    scalar recTimeStep_=recurrenceBase.recM().recTimeStep();
+
+    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,"Flow");
+        #include "TEq.H"
+        particleCloud.clockM().stop("Flow");
+
+
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+            Info << "Updating fields at run time " << runTime.timeOutputValue()
+                 << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            recurrenceBase.updateRecFields();
+            #include "readFields.H"
+            recTimeIndex++;
+        }
+
+        particleCloud.clockM().start(27,"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/testTwoFluidRecurrenceTurbulence/Allwclean

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wclean libso recurrenceTurbulence
+
+wclean
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/testTwoFluidRecurrenceTurbulence/Allwmake

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wmake libso recurrenceTurbulence
+
+wmake
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/testTwoFluidRecurrenceTurbulence/CEqn.H

@@ -0,0 +1,29 @@
+// build equation system
+/*
+    Note the use of the effective viscosity, which is provided by the turbulence model
+    The recurrence-based turbulence models are derived from the standard base classes 
+    of OpenFOAM, thus they behave as a normal turbulence model would.
+*/
+
+alphaRhoPhiCarrier = linearInterpolate(alpha2*rhoCarrier)*phi2;
+
+fvScalarMatrix CEqn
+(
+	fvm::ddt(alphaCarrier*rhoCarrier, C)
+  + fvm::div(alphaRhoPhiCarrier, C, "div(alphaRhoPhi,C)")
+  - fvm::Sp(fvc::div(alphaRhoPhiCarrier), C)
+
+  - fvm::laplacian
+	(
+		fvc::interpolate(alpha2)
+	   *fvc::interpolate(carrierPhase.turbulence().muEff()/Sc),
+		C
+	)
+	==
+	fvm::SuSp(alphaCarrier*(1.0 - alphaCarrier)*rhoCarrier*K, C)
+	+ fvOptions(alphaCarrier*rhoCarrier, C)
+);
+
+// solve equations
+fvOptions.constrain(CEqn);
+CEqn.solve();

applications/solvers/testTwoFluidRecurrenceTurbulence/Make/files

@@ -0,0 +1,3 @@
+testTwoFluidRecurrenceTurbulence.C
+
+EXE = $(FOAM_USER_APPBIN)/testTwoFluidRecurrenceTurbulence

applications/solvers/testTwoFluidRecurrenceTurbulence/Make/options

@@ -0,0 +1,31 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/phaseSystems/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/interfacialModels/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/interfacialCompositionModels/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/twoPhaseCompressibleTurbulenceModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I../../../src/recurrence/lnInclude \
+    -IrecurrenceTurbulence/lnInclude
+
+EXE_LIBS = \
+    -lreactingPhaseSystem \
+    -lreactingTwoPhaseSystem \
+    -lreactingEulerianInterfacialModels \
+    -lreactingEulerianInterfacialCompositionModels \
+    -ltwoPhaseReactingTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -L$(FOAM_USER_LIBBIN) \
+    -lrecurrence \
+    -lrecurrenceTwoPhaseTurbulenceModels

applications/solvers/testTwoFluidRecurrenceTurbulence/checkTurbulenceModels.H

@@ -0,0 +1,70 @@
+//===============================
+// recurrence turbulence
+//===============================
+
+
+    // check both phases for turbulence models
+    forAllIter(PtrListDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        
+        Info << "Checking phase " << phase.name() << "'s turbulence model: "
+            << phase.turbulence().type() << endl;
+        
+        /*
+            Check for laminar turbulence. This works with OpenFOAM-4.0 and OpenFOAM-5.0,
+            as the laminar, multi-phase turbulence model is named "laminar" in OF-4.0
+            and "Stokes" in OF-5.0
+        */
+        if (phase.turbulence().type() == "laminar" || phase.turbulence().type() == "Stokes")
+        {
+            // do nothing
+        }
+        else if (isA<Foam::recurrenceTurbulenceModel>(phase.turbulence()))
+        {
+            /*
+                create a reference of the type recurrenceTurbulenceModel
+                register the recurrence model with the recurrenceTurbulenceModel
+            */
+            
+            // get const-reference to the turbulence model
+            const phaseCompressibleTurbulenceModel& turbConstRef = phase.turbulence();
+            
+            // cast away const-ness, the underlying turbulence model is not a const object, so this is bad but fine
+            phaseCompressibleTurbulenceModel& turbRef = const_cast<phaseCompressibleTurbulenceModel&>(turbConstRef);
+            
+            // cast away the wrapper class, to get a reference to the turbulence models' base class
+            PhaseCompressibleTurbulenceModel<phaseModel>& baseTurbRef
+            (
+                static_cast<PhaseCompressibleTurbulenceModel<phaseModel>&>(turbRef)
+            );
+            
+            // casting down the family tree
+            Foam::recurrenceTurbulenceModel& recTurbRef
+            (
+                dynamic_cast<Foam::recurrenceTurbulenceModel&>(baseTurbRef)
+            );
+            
+            
+            // set recurrenceBase pointer
+            recTurbRef.setRecurrenceBasePtr(&recurrenceBase);
+            
+            
+            // check model settings
+            turbRef.validate();
+        }
+        else
+        {
+            /*
+                In a recurrence run, we do not compute any turbulence as we do not solve the fluid flow
+                At this point, the phase is not laminar (i.e. not using turbulence) or 
+                    using recurrenceTurbulence (i.e. taking turbulent quantities from the data base).
+                Hence, abort!
+            */
+            FatalError
+                << "Wrong turbulence model type "
+                << phase.turbulence().type() << " for phase " << phase.name() << nl << nl
+                << "Valid turbulence model types are types derived from recurrenceTurbulenceModel or laminar" << endl
+                << exit(FatalError);
+        }
+    }

applications/solvers/testTwoFluidRecurrenceTurbulence/createTransportFields.H

@@ -0,0 +1,79 @@
+/* --------------------------------------------------------------------------------- */
+/* 								read flotation properties							 */
+/* --------------------------------------------------------------------------------- */
+
+Info<< "Reading scalarTransportProperties\n" << endl;
+IOdictionary scalarTransportProperties
+(
+    IOobject
+    (
+        "scalarTransportProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+
+const scalar Sc(scalarTransportProperties.lookupOrDefault<scalar>("Sc",scalar(1.0)));
+
+
+const word carrierPhaseName(scalarTransportProperties.lookup("carrierPhase"));
+
+if (carrierPhaseName != phase1.name() && carrierPhaseName != phase2.name())
+{
+	FatalError << "No valid carrier phase specified" << nl
+		<< "Valid phase names are: " << nl
+		<< phase1.name() << ",  " << phase2.name()
+		<< abort(FatalError);
+}
+
+phaseModel& carrierPhase = (carrierPhaseName == phase1.name()) ? phase1 : phase2;
+
+const word dispersePhaseName = (carrierPhaseName == phase1.name()) ? phase2.name() : phase1.name();
+
+
+volScalarField& rhoCarrier = carrierPhase.thermo().rho();
+volScalarField& alphaCarrier = carrierPhase;
+surfaceScalarField& alphaRhoPhiCarrier = carrierPhase.alphaRhoPhi();
+
+
+volScalarField contErrCarrier
+(
+    "contErrCarrier",
+    fvc::ddt(alphaCarrier, rhoCarrier)
+);
+
+
+
+Info<< "Reading field C\n" << endl;
+volScalarField C
+(
+    IOobject
+    (
+        "C",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+volScalarField K
+(
+    IOobject
+    (
+        "K",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+
+
+

applications/solvers/testTwoFluidRecurrenceTurbulence/readFields.H

@@ -0,0 +1,9 @@
+// update flow fields
+recurrenceBase.recM().exportVolScalarField("alpha."+carrierPhaseName,alpha2);
+recurrenceBase.recM().exportVolScalarField("alpha."+dispersePhaseName,alpha1);
+
+recurrenceBase.recM().exportVolVectorField("U."+carrierPhaseName,U2);
+
+// update turbulence models
+phase1.correctTurbulence();
+phase2.correctTurbulence();

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/Allwclean

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

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/Allwmake

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

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/Make/files

@@ -0,0 +1,10 @@
+recurrenceTurbulenceModel/recurrenceTurbulenceModel.C
+
+
+recurrenceTurbulenceModels.C
+
+recurrenceKEpsilon/recurrenceKEpsilon.C
+recurrenceKOmega/recurrenceKOmega.C
+recurrenceSmagorinsky/recurrenceSmagorinsky.C
+
+LIB = $(FOAM_USER_LIBBIN)/librecurrenceTwoPhaseTurbulenceModels

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/Make/options

@@ -0,0 +1,27 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/phaseSystems/lnInclude \
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/interfacialModels/lnInclude\
+    -I$(FOAM_SOLVERS)/multiphase/reactingEulerFoam/interfacialCompositionModels/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/transportModel \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I../recurrenceTurbulenceModel/lnInclude \
+    -I../../../../src/recurrence/lnInclude
+
+LIB_LIBS = \
+    -lreactingPhaseSystem \
+    -lreactingTwoPhaseSystem \
+    -lreactingEulerianInterfacialModels \
+    -lreactingEulerianInterfacialCompositionModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -L$(FOAM_USER_LIBBIN) \
+    -lrecurrence \
+    -lreactingTwoPhaseSystem

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/ReadMe.md

@@ -0,0 +1,107 @@
+# Recurrence-based, multi-phase turbulence modelling
+
+This model implements recurrence-based turbulence models, i.e. the fundamental 
+turbulent field quantities are read from the data base and are not solved for.
+All derived field quantities are computed just in the same way as the proper 
+turbulence models do. By deriving the recurrence-based turbulence models from 
+somewhere up the family tree of OpenFOAM's turbulence model class hierarchy, 
+the recurrence-based turbulence models are fully compatible with OpenFOAM's 
+generic treatment of turbulence modelling, i.e. solvers and libraries interact 
+with references to a generic base type of the actual turbulence model. Hence, 
+solvers and libraries may remain blissfully ignorant of the actual turbulence 
+model in use.
+
+For laminar phases no special treatment is necessary, as the *laminar* 
+turbulence model does not compute any fields.
+
+
+## Development notes
+
+The initial development covers only a small number of turbulence models.
+
+
+## Notes on usage
+
+The turbulence model in use for the recurrence run must be the recurrence-based 
+equivalent of the turbulence model used for generating the data base, i.e. if 
+the data base was computed using the *kEpsilon* model, then the recurrence solver 
+is to employ the *recurrenceKEpsilon* turbulence model. This model will read 
+the relevant model coefficients from the *turbulenceProperties* dictionary, and 
+make sure that the turbulent fields `k` and `epsilon` are contained in the data 
+base.
+
+Whenever, the solver or a library calls `turbulence->nut()` to access the 
+turbulent viscosity, the recurrence-based kEpsilon model will compute `nut` 
+according to kEpsilon's relations `nut = Cmu*sqr(k)/epsilon`, with the fields 
+`k` and `epsilon` being from the current snapshot provided by the recurrence model.
+
+Thus, the fundamental turbulent field quantities of the employed turbulence model 
+have to be added to the *volScalarFields* list in the `recProperties` dictionary 
+controlling the recurrence model. This will ensure that the turbulent field 
+quantities are read from the data base.
+
+
+## Notes on the implementation
+
+The base class implements the method `void setRecurrenceBasePtr(recBase*)`, which 
+is used to give the recurrence-based turbulence models a reference (technically 
+a pointer) to the recurrence model. Thus, after construction of the turbulence 
+models and the recurrence model, `setRecurrenceBasePtr()` needs to be called as 
+the pointer to the recurrence model is initialized by the constructor with `NULL`.
+Trying to access the recurrence model from within the recurrence-based turbulence 
+model prior to setting the pointer to the recurrence model with 
+`setRecurrenceBasePtr()` will result in a segmentation fault.
+In order to be able to call `setRecurrenceBasePtr()`, the generic reference to 
+the turbulence model needs to be converted into a reference of the base class' 
+type, i.e. `recurrenceTurbulenceModel`.
+This unfortunate deviation from good standards, i.e. making full use of C++'s 
+polymorphism, should be the only instance of having to use non-pretty hacks.
+However, apart from initialisation, i.e. setting the pointer to the recurrence 
+model, the recurrence-based turbulence models adhere to the generic interface of 
+OpenFOAM's turbulence models, and can be used as any other turbulence model.
+
+
+The concrete implementations, e.g. *recurrenceKEpsilon*, use the method 
+`validate()` to check whether the underlying turbulent quantities are specified 
+for use in the data base in the *volScalarFields* list in the `recProperties` 
+dictionary. This method is part of the signature of the class `Foam::turbulenceModel`, 
+which is the very base class of all turbulence models in OpenFOAM.
+In proper turbulence models, this method is used to check whether the internal 
+fields are properly initialized and to update all derived quantities.
+In the solver, `validate()` must not be called prior to `setRecurrenceBasePtr()`, 
+as validate accesses the recurrence model. The wrong order of function calls will 
+result in a segmentation fault, as the pointer to the recurrence model is 
+initialized by the constructor with `NULL`.
+
+
+The concrete implementations, e.g. *recurrenceKEpsilon*, use the method 
+`correct()` to update the turbulent field quantities from the data base, 
+and in turn update the derived quantities, such as `nut`.
+This method is part of the signature of the class `Foam::turbulenceModel`, 
+which is the very base class of all turbulence models in OpenFOAM.
+In proper turbulence models, this method is used to solve for the next time step.
+
+
+## Compilation
+
+Source OpenFOAM and simply compile with
+
+```bash
+./Allwclean
+./Allwmake
+```
+
+The script `Allwclean` will clear all previous builds. This step is not needed for
+first-time compilation. It is, however, recommended for subsequent compilations, as
+it completely clears the slate. The script `Allwmake` will run the compilation for
+the passive particle model.
+
+
+## Required software
+
+This model has been tested with the following versions of OpenFOAM:
+
+* OpenFOAM-4.0
+* OpenFOAM-5.0
+
+

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceKEpsilon/recurrenceKEpsilon.C

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "recurrenceKEpsilon.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+void Foam::RASModels::recurrenceKEpsilon::correctNut()
+{
+    this->nut_ = Cmu_*sqr(k_)/epsilon_;
+    this->nut_.correctBoundaryConditions();
+}
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::RASModels::recurrenceKEpsilon::recurrenceKEpsilon
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& phase,
+    const word& propertiesName,
+    const word& type
+
+)
+:
+    eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        phase,
+        propertiesName
+    ),
+    recurrenceTurbulenceModel(U.group()),
+    Cmu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar("k0", dimensionSet(0,2,-2,0,0), 0.0)
+    ),
+    epsilon_
+    (
+        IOobject
+        (
+            IOobject::groupName("epsilon", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar("eps0", dimensionSet(0,2,-3,0,0), 0.0)
+    )
+{
+    if (type == typeName)
+    {
+        printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+Foam::RASModels::recurrenceKEpsilon::~recurrenceKEpsilon()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::RASModels::recurrenceKEpsilon::read()
+{
+    if
+    (
+        eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>::read()
+    )
+    {
+        Cmu_.readIfPresent(this->coeffDict());
+        
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::RASModels::recurrenceKEpsilon::correct()
+{
+    // update turbulence fields
+    recurrenceBasePtr_->recM().exportVolScalarField("k."+group_, this->k_);
+    recurrenceBasePtr_->recM().exportVolScalarField("epsilon."+group_, this->epsilon_);
+    
+    // update nut
+    correctNut();
+}
+
+
+void Foam::RASModels::recurrenceKEpsilon::validate()
+{
+    /*
+        Check whether k and epsilon are included in the dataBase.
+        The check only makes sure that these fields are included in the
+            volScalarFields list of recProperties.
+        Whether the fields are actually contained in the dataBase is 
+            done by the recurrenceModel itself.
+    */
+    bool foundK(false);
+    bool foundEpsilon(false);
+    
+    wordList fieldNames(recurrenceBasePtr_->recM().volScalarFieldNames());
+    
+    forAll(fieldNames, i)
+    {
+        word curFieldName = fieldNames[i];
+        
+        if (curFieldName == k_.name())
+        {
+            foundK = true;
+        }
+        
+        if (curFieldName == epsilon_.name())
+        {
+            foundEpsilon = true;
+        }
+        
+        
+    }
+    
+    if (not (foundK and foundEpsilon))
+    {
+        FatalError
+            << "Fields " << k_.name() << " and " << epsilon_.name()
+            << " not specified in the volScalarFields list of recProperties!" << nl
+            << "volScalarFields : " << fieldNames << nl
+            << "Add these fields and make sure they are contained in the dataBase." << nl
+            << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceKEpsilon/recurrenceKEpsilon.H

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::recurrenceKEpsilon
+
+Description
+    recurrence-based kEpsilon turbulence model
+
+    This model provides kEpslion's turbulence quantities that were computed 
+    elsewhere, i.e. taken from the recurrence dataBase.
+    To be used by recurrence solvers.
+
+SourceFiles
+    recurrenceKEpsilon.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef recurrenceKEpsilon_H
+#define recurrenceKEpsilon_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "EddyDiffusivity.H"
+
+#include "recurrenceTurbulenceModel.H"
+
+#include "autoPtr.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class recurrenceKEpsilon Declaration
+\*---------------------------------------------------------------------------*/
+
+class recurrenceKEpsilon
+:
+    public eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>,
+    public recurrenceTurbulenceModel
+{
+    // Private data
+
+    // Private Member Functions
+
+
+        //- Disallow default bitwise copy construct
+        recurrenceKEpsilon(const recurrenceKEpsilon&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const recurrenceKEpsilon&);
+        
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar Cmu_;
+
+        // Fields
+
+            volScalarField k_;
+            volScalarField epsilon_;
+    
+    // Protected Member Functions
+
+        virtual void correctNut();
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("recurrenceKEpsilon");
+
+
+    // Constructors
+
+        //- Construct from components
+        recurrenceKEpsilon
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const phaseModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+
+    //- Destructor
+    virtual ~recurrenceKEpsilon();
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+        
+        
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return epsilon_;
+        }
+
+        //- Update the turbulent fields
+        virtual void correct();
+        
+        //- Check model settings
+        virtual void validate();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceKOmega/recurrenceKOmega.C

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "recurrenceKOmega.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+void Foam::RASModels::recurrenceKOmega::correctNut()
+{
+    this->nut_ = k_/omega_;
+    this->nut_.correctBoundaryConditions();
+}
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::RASModels::recurrenceKOmega::recurrenceKOmega
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& phase,
+    const word& propertiesName,
+    const word& type
+
+)
+:
+    eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        phase,
+        propertiesName
+    ),
+    recurrenceTurbulenceModel(U.group()),
+    Cmu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar("k0", dimensionSet(0,2,-2,0,0), 0.0)
+    ),
+    omega_
+    (
+        IOobject
+        (
+            IOobject::groupName("omega", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar("om0", dimensionSet(0,0,-1,0,0), 0.0)
+    )
+{
+    if (type == typeName)
+    {
+        printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+Foam::RASModels::recurrenceKOmega::~recurrenceKOmega()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::RASModels::recurrenceKOmega::read()
+{
+    if
+    (
+        eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>::read()
+    )
+    {
+        Cmu_.readIfPresent(this->coeffDict());
+        
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::RASModels::recurrenceKOmega::correct()
+{
+    // update turbulence fields
+    recurrenceBasePtr_->recM().exportVolScalarField("k."+group_, this->k_);
+    recurrenceBasePtr_->recM().exportVolScalarField("omega."+group_, this->omega_);
+    
+    // update nut
+    correctNut();
+}
+
+
+void Foam::RASModels::recurrenceKOmega::validate()
+{
+    /*
+        Check whether k and omega are included in the dataBase.
+        The check only makes sure that these fields are included in the
+            volScalarFields list of recProperties.
+        Whether the fields are actually contained in the dataBase is 
+            done by the recurrenceModel itself.
+    */
+    bool foundK(false);
+    bool foundOmega(false);
+    
+    wordList fieldNames(recurrenceBasePtr_->recM().volScalarFieldNames());
+    
+    forAll(fieldNames, i)
+    {
+        word curFieldName = fieldNames[i];
+        
+        if (curFieldName == k_.name())
+        {
+            Info << "Found " << k_.name()<< endl;
+            foundK = true;
+        }
+        
+        if (curFieldName == omega_.name())
+        {
+            Info << "Found " << omega_.name()<< endl;
+            foundOmega = true;
+        }
+    }
+    
+    if (not (foundK and foundOmega))
+    {
+        FatalError
+            << "Fields " << k_.name() << " and " << omega_.name()
+            << " not specified in the volScalarFields list of recProperties!" << nl
+            << "volScalarFields : " << fieldNames << nl
+            << "Add these fields and make sure they are contained in the dataBase." << nl
+            << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceKOmega/recurrenceKOmega.H

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::recurrenceKOmega
+
+Description
+    recurrence-based kOmega turbulence model
+
+    This model provides kOmega's turbulence quantities that were computed 
+    elsewhere, i.e. taken from the recurrence dataBase
+    To be used by recurrence solvers
+
+SourceFiles
+    recurrenceKOmega.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef recurrenceKOmega_H
+#define recurrenceKOmega_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "EddyDiffusivity.H"
+
+#include "recurrenceTurbulenceModel.H"
+
+#include "autoPtr.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class recurrenceKOmega Declaration
+\*---------------------------------------------------------------------------*/
+
+class recurrenceKOmega
+:
+    public eddyViscosity<RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>>,
+    public recurrenceTurbulenceModel
+{
+    // Private data
+
+    // Private Member Functions
+
+
+        //- Disallow default bitwise copy construct
+        recurrenceKOmega(const recurrenceKOmega&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const recurrenceKOmega&);
+        
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar Cmu_;
+
+        // Fields
+
+            volScalarField k_;
+            volScalarField omega_;
+    
+    // Protected Member Functions
+
+        virtual void correctNut();
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("recurrenceKOmega");
+
+
+    // Constructors
+
+        //- Construct from components
+        recurrenceKOmega
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const phaseModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+
+    //- Destructor
+    virtual ~recurrenceKOmega();
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+        
+        
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulence specific dissipation rate
+        virtual tmp<volScalarField> omega() const
+        {
+            return omega_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "epsilon",
+                        this->mesh_.time().timeName(),
+                        this->mesh_
+                    ),
+                    Cmu_*k_*omega_,
+                    omega_.boundaryField().types()
+                )
+            );
+        }
+
+        //- Update the turbulent fields
+        virtual void correct();
+        
+        //- Check model settings
+        virtual void validate();
+        
+        
+    // Setters
+        
+        void setRecurrenceBasePtr(recBase* recurrenceBasePtr)
+        {
+            recurrenceBasePtr_ = recurrenceBasePtr;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceSmagorinsky/recurrenceSmagorinsky.C

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "recurrenceSmagorinsky.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::LESModels::recurrenceSmagorinsky::k
+(
+    const tmp<volTensorField>& gradU
+) const
+{
+    volSymmTensorField D(symm(gradU));
+
+    volScalarField a(this->Ce_/this->delta());
+    volScalarField b((2.0/3.0)*tr(D));
+    volScalarField c(2*Ck_*this->delta()*(dev(D) && D));
+
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("k", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_
+            ),
+            sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a))
+        )
+    );
+}
+
+
+void Foam::LESModels::recurrenceSmagorinsky::correctNut()
+{
+    volScalarField k(this->k(fvc::grad(this->U_)));
+
+    this->nut_ = Ck_*this->delta()*sqrt(k);
+    this->nut_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::LESModels::recurrenceSmagorinsky::recurrenceSmagorinsky
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& phase,
+    const word& propertiesName,
+    const word& type
+)
+:
+    LESeddyViscosity<EddyDiffusivity<phaseCompressibleTurbulenceModel>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        phase,
+        propertiesName
+    ),
+    recurrenceTurbulenceModel(U.group()),
+    Ck_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ck",
+            this->coeffDict_,
+            0.094
+        )
+    )
+{
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::LESModels::recurrenceSmagorinsky::read()
+{
+    if (LESeddyViscosity<EddyDiffusivity<phaseCompressibleTurbulenceModel>>::read())
+    {
+        Ck_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::LESModels::recurrenceSmagorinsky::epsilon() const
+{
+    volScalarField k(this->k(fvc::grad(this->U_)));
+
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("epsilon", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_
+            ),
+            this->Ce_*k*sqrt(k)/this->delta()
+        )
+    );
+}
+
+
+void Foam::LESModels::recurrenceSmagorinsky::correct()
+{
+    // update nut
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceSmagorinsky/recurrenceSmagorinsky.H

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::recurrenceSmagorinsky
+
+Description
+    The recurrenceSmagorinsky SGS model.
+
+    This model provides Smagorinsky's turbulence quantities.
+    To be used by recurrence solvers.
+
+SourceFiles
+    recurrenceSmagorinsky.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef recurrenceSmagorinsky_H
+#define recurrenceSmagorinsky_H
+
+#include "LESModel.H"
+#include "LESeddyViscosity.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "EddyDiffusivity.H"
+
+#include "recurrenceTurbulenceModel.H"
+
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class recurrenceSmagorinsky Declaration
+\*---------------------------------------------------------------------------*/
+
+class recurrenceSmagorinsky
+:
+    public LESeddyViscosity<EddyDiffusivity<phaseCompressibleTurbulenceModel>>,
+    public recurrenceTurbulenceModel
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        recurrenceSmagorinsky(const recurrenceSmagorinsky&);
+        void operator=(const recurrenceSmagorinsky&);
+
+
+protected:
+
+    // Protected data
+
+        dimensionedScalar Ck_;
+
+
+    // Protected Member Functions
+
+        //- Return SGS kinetic energy
+        //  calculated from the given velocity gradient
+        tmp<volScalarField> k(const tmp<volTensorField>& gradU) const;
+
+        //- Update the SGS eddy viscosity
+        virtual void correctNut();
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("recurrenceSmagorinsky");
+
+
+    // Constructors
+
+        //- Construct from components
+        recurrenceSmagorinsky
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const phaseModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~recurrenceSmagorinsky()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return SGS kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k(fvc::grad(this->U_));
+        }
+
+        //- Return sub-grid disipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Correct Eddy-Viscosity and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceTurbulenceModel/recurrenceTurbulenceModel.C

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "recurrenceTurbulenceModel.H"
+
+
+namespace Foam
+{
+    defineTypeNameAndDebug(recurrenceTurbulenceModel, 0);
+}
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::recurrenceTurbulenceModel::recurrenceTurbulenceModel
+(
+    const word group
+)
+:
+    recurrenceBasePtr_(NULL),
+    group_(group)
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::recurrenceTurbulenceModel::~recurrenceTurbulenceModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceTurbulenceModel/recurrenceTurbulenceModel.H

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::recurrenceTurbulenceModel
+
+Description
+    recurrence-based turbulence model base class
+
+    This turbulence model class provides a framework for derived turbulence
+    models to provide turbulence quantities that were computed elsewhere, 
+    i.e. taken from the recurrence dataBase.
+    
+    A concrete recurrence-based turbulence model implementation, e.g. recurrenceKEpsilon, 
+    needs to be derived from the base class of the original turbulence model 
+    and this class. This class provides the link to the recurrence model. The original 
+    base class provides the proper turbulence modelling interfaces of OpenFOAM's 
+    turbulence modelling framework.
+    
+    To be used by recurrence solvers.
+
+SourceFiles
+    recurrenceTurbulenceModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef recurrenceTurbulenceModel_H
+#define recurrenceTurbulenceModel_H
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "autoPtr.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class recurrenceTurbulenceModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class recurrenceTurbulenceModel
+{
+    // Private Member Functions
+
+
+        //- Disallow default bitwise copy construct
+        recurrenceTurbulenceModel(const recurrenceTurbulenceModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const recurrenceTurbulenceModel&);
+        
+
+protected:
+
+    // Protected data
+
+        recBase* recurrenceBasePtr_;
+        const word group_;
+        
+    
+    // Protected Member Functions
+
+
+public:
+
+    //- Runtime type information
+    TypeName("recurrenceTurbulenceModel");
+
+
+    // Constructors
+
+        //- Construct from components
+        recurrenceTurbulenceModel
+        (
+            const word group
+        );
+
+
+    //- Destructor
+    virtual ~recurrenceTurbulenceModel();
+
+
+    // Member Functions
+        
+        void setRecurrenceBasePtr(recBase* recurrenceBasePtr)
+        {
+            recurrenceBasePtr_ = recurrenceBasePtr;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/recurrenceTurbulence/recurrenceTurbulenceModels.C

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phaseCompressibleTurbulenceModel.H"
+#include "addToRunTimeSelectionTable.H"
+#include "makeTurbulenceModel.H"
+
+#include "laminar.H"
+#include "RASModel.H"
+#include "LESModel.H"
+
+#include "recurrenceKEpsilon.H"
+#include "recurrenceKOmega.H"
+#include "recurrenceSmagorinsky.H"
+
+
+
+// Instructions for OpenFOAM-5.0 
+/*makeTurbulenceModelTypes
+(
+    volScalarField,
+    volScalarField,
+    compressibleTurbulenceModel,
+    PhaseCompressibleTurbulenceModel,
+    ThermalDiffusivity,
+    phaseModel
+);
+
+
+makeTurbulenceModel
+(phaseModelPhaseCompressibleTurbulenceModel, RAS, recurrenceKEpsilon);
+
+makeTurbulenceModel
+(phaseModelPhaseCompressibleTurbulenceModel, RAS, recurrenceKOmega);
+
+makeTurbulenceModel
+(phaseModelPhaseCompressibleTurbulenceModel, LES, recurrenceSmagorinsky);
+*/
+
+
+
+// Instructions for OpenFOAM-4.0 
+makeBaseTurbulenceModel
+(
+    volScalarField,
+    volScalarField,
+    compressibleTurbulenceModel,
+    PhaseCompressibleTurbulenceModel,
+    ThermalDiffusivity,
+    phaseModel
+);
+
+#define makeRASModel(Type)                                            \
+    makeTurbulenceModel                                               \
+    (phaseModelPhaseCompressibleTurbulenceModel, RAS, Type)
+
+#define makeLESModel(Type)                                            \
+    makeTurbulenceModel                                               \
+    (phaseModelPhaseCompressibleTurbulenceModel, LES, Type)
+
+makeRASModel(recurrenceKEpsilon);
+
+makeRASModel(recurrenceKOmega);
+
+makeLESModel(recurrenceSmagorinsky);
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/testTwoFluidRecurrenceTurbulence.C

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    testTwoFluidRecurrenceTurbulence
+
+Description
+    A modified variant of the two-fluid, recurrence model A solver
+    with the extension of recurrence-based, multi-phase turbulence modelling.
+    This application is used to test whether turbulent fields can be provided
+    by the recurrence-based turbulence models.
+    
+    Run this test application in a recurrence case, with turbulence enabled and 
+    the necessary turbulent field quantities present in the data base.
+    
+    Note the initialisation in checkTurbulenceModels.H
+    Updating the turbulence model is done by calling phaseX.correctTurbulence()
+    in the file readFields.H
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "twoPhaseSystem.H"
+#include "phaseCompressibleTurbulenceModel.H"
+#include "pimpleControl.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "recurrenceTurbulenceModel.H"
+
+/* // uncomment for OpenFOAM-5.0
+namespace Foam
+{
+    tmp<volScalarField> byDt(const volScalarField& vf)
+    {
+        if (fv::localEulerDdt::enabled(vf.mesh()))
+        {
+            return fv::localEulerDdt::localRDeltaT(vf.mesh())*vf;
+        }
+        else
+        {
+            return vf/vf.mesh().time().deltaT();
+        }
+    }
+
+    tmp<surfaceScalarField> byDt(const surfaceScalarField& sf)
+    {
+        if (fv::localEulerDdt::enabled(sf.mesh()))
+        {
+            return fv::localEulerDdt::localRDeltaTf(sf.mesh())*sf;
+        }
+        else
+        {
+            return sf/sf.mesh().time().deltaT();
+        }
+    }
+}
+*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H" // remove for OpenFOAM-5.0
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    
+    #include "createTransportFields.H"
+
+    if (!LTS)
+    {
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+    
+    Switch faceMomentum
+    (
+        pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
+    );
+
+    
+    recBase recurrenceBase(mesh);
+    
+    #include "checkTurbulenceModels.H"
+
+    #include "pUf/createDDtU.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting recurrence-based time loop\n" << endl;
+    
+    label recTimeIndex(0);
+    scalar recTimeStep_=recurrenceBase.recM().recTimeStep();
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        #include "CourantNos.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        
+        #include "CEqn.H"
+        
+        
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+	        Info << "Updating fields at run time " << runTime.timeOutputValue()
+	            << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+	        
+            recurrenceBase.updateRecFields();
+	        #include "readFields.H"
+            recTimeIndex++;
+        }
+        
+
+        runTime.write();
+        
+        #include "writeCField.H"
+
+        Info<< "ExecutionTime = "
+            << runTime.elapsedCpuTime()
+            << " s\n\n" << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/testTwoFluidRecurrenceTurbulence/writeCField.H

@@ -0,0 +1,12 @@
+/* ----------------------------------------------
+	Write averaged particle volume fraction
+---------------------------------------------- */
+
+// essential information
+Info << "Total mass :";
+Info << tab << sum(C*rhoCarrier*alphaCarrier*mesh.V());
+Info << endl;
+
+Info << "Total Carrier mass :";
+Info << tab << sum(rhoCarrier*alphaCarrier*mesh.V());
+Info << endl;

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/mapDataBase/mapFieldData.sh

@@ -0,0 +1,30 @@
+#!/bin/sh
+# Source run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# to be executed from one level above the source directory
+
+if [ $# -eq 0 ]
+  then
+    sourceName="dataBase"
+    targetName="dataBaseCoarse"
+  else
+    sourceName=$1
+    targetName=$2
+fi
+
+cd $sourceName
+
+for time in *
+do
+    if [ $time != "system" ] && [ $time != "constant" ];
+    then
+        cd ../$targetName
+        echo "Found $time."
+        sed -i "/^startTime/c\startTime \t$time;" ./system/controlDict
+        grep 'startTime' ./system/controlDict
+        mapFields ../$sourceName -sourceTime $time -consistent
+        cd ../$sourceName
+    fi
+done
+

applications/utilities/rBaseMirror/mirrorProperties

@@ -0,0 +1,39 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          mirrorProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//===========================================================================//
+// sub-models & settings
+
+
+refPoint (0 0 0);
+
+refDirection (1 0 0);
+
+fieldName U;
+
+dataBaseName dataBase;
+
+
+// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirrorScalar/Make/files

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

applications/utilities/rBaseMirror/rBaseMirrorScalar/Make/options

@@ -0,0 +1,17 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions

applications/utilities/rBaseMirror/rBaseMirrorScalar/createFields.H

@@ -0,0 +1,17 @@
+    IOdictionary mirrorProperties
+    (
+        IOobject
+        (
+            "mirrorProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    vector refPoint(mirrorProperties.lookup("refPoint"));
+    vector refDirection(mirrorProperties.lookup("refDirection"));
+
+    word fieldName(mirrorProperties.lookup("fieldName"));
+    word dataBaseName(mirrorProperties.lookup("dataBaseName"));

applications/utilities/rBaseMirror/rBaseMirrorScalar/rBaseMirrorScalar.C

@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rBaseMirror
+
+
+Description
+    Read time series and extend it by mirrored fields if geometry possesses
+    the same symmetry
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    timeSelector::addOptions();
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+   // read in start and end time from controlDict
+
+    scalar startTime=runTime.startTime().value();
+    scalar endTime=runTime.endTime().value();
+    scalar origTimeRange = endTime - startTime;
+
+    Info << "start time = " << runTime.startTime() << endl;
+    Info << "end time = " << runTime.endTime() << endl;
+
+    // check which time directories are present
+    // instantList timeDirs = timeSelector::select0(runTime, args);
+    // runTime.setTime(timeDirs[0], 0);
+
+    #include "createMesh.H"
+
+    #include "createFields.H"
+
+    Foam::Time recTime(fileName(dataBaseName), "", "../system", "../constant", false);
+    instantList timeDirs(recTime.times());
+    recTime.setTime(timeDirs[0],0);
+
+    #include "readFields.H"
+
+    Info << fieldName << endl;
+
+    volScalarField transformedField = origField;
+
+    scalar t;
+
+    label shiftedTimeI = 0;
+
+    // check number of time directories
+    label shift = 0;
+    forAll(timeDirs, timeI)
+    {
+
+       if (recTime.timeName() == "constant") continue;
+       recTime.setTime(timeDirs[timeI], timeI);
+       t = recTime.value();
+       if(t < startTime) continue;
+       if(t > endTime) continue;
+       shift++;
+    }
+
+    scalar dt = origTimeRange / (shift - 1.0);
+    recTime.setEndTime(startTime + 2 * origTimeRange + dt);
+
+    label cellI_transformed = -1;
+    forAll(timeDirs, timeI)
+    {
+
+       recTime.setTime(timeDirs[timeI], timeI);
+       t = recTime.value();
+       if(t < startTime) continue;
+       if(t > endTime) continue;
+       Info << "time = " << t << ", time index = " << timeI << endl;
+
+       #include "readFields.H"
+
+       forAll(transformedField, cellI)
+       {
+           vector position = mesh.C()[cellI];
+           vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
+           cellI_transformed = mesh.findCell(transformedPosition);
+           if(cellI_transformed < 0)
+           {
+               Info << "Couldn't find transformed cell. Stopping." << endl;
+               return 0;
+           }
+
+           scalar value = origField[cellI_transformed];
+           scalar transformedValue = value;
+
+           transformedField[cellI] = transformedValue;
+       }
+
+       shiftedTimeI = timeI + shift;
+       t = recTime.value() + origTimeRange + dt;
+       runTime.setTime(t, shiftedTimeI);
+       Info << "creating transformed fields for time = " << t << ", time index = " << shiftedTimeI << endl;
+       transformedField.write();
+    }
+
+    Info << "\nEnd" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirrorScalar/readFields.H

@@ -0,0 +1,13 @@
+
+    volScalarField origField
+    (
+        IOobject
+        (
+            fieldName,
+            recTime.timePath(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh
+    );

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/files

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

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/options

@@ -0,0 +1,17 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions

applications/utilities/rBaseMirror/rBaseMirrorVec/createFields.H

@@ -0,0 +1,17 @@
+    IOdictionary mirrorProperties
+    (
+        IOobject
+        (
+            "mirrorProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    vector refPoint(mirrorProperties.lookup("refPoint"));
+    vector refDirection(mirrorProperties.lookup("refDirection"));
+
+    word fieldName(mirrorProperties.lookup("fieldName"));
+    word dataBaseName(mirrorProperties.lookup("dataBaseName"));

applications/utilities/rBaseMirror/rBaseMirrorVec/rBaseMirrorVec.C

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rBaseMirror
+
+Description
+    Read time series and extend it by mirrored fields if geometry possesses
+    the same symmetry
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    timeSelector::addOptions();
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+   // read in start and end time from controlDict
+
+    scalar startTime=runTime.startTime().value();
+    scalar endTime=runTime.endTime().value();
+    scalar origTimeRange = endTime - startTime;
+
+    Info << "start time = " << runTime.startTime() << endl;
+    Info << "end time = " << runTime.endTime() << endl;
+
+    // check which time directories are present
+    //instantList timeDirs = timeSelector::select0(runTime, args);
+    //runTime.setTime(timeDirs[0], 0);
+
+    #include "createMesh.H"
+
+    #include "createFields.H"
+
+    Foam::Time recTime(fileName(dataBaseName), "", "../system", "../constant", false);
+    instantList timeDirs(recTime.times());
+    recTime.setTime(timeDirs[0],0);
+
+    #include "readFields.H"
+
+    Info << fieldName << endl;
+
+    volVectorField transformedField = origField;
+
+    scalar t;
+
+    label shiftedTimeI = 0;
+
+    // check number of time directories
+    label shift = 0;
+    forAll(timeDirs, timeI)
+    {
+
+       if (recTime.timeName() == "constant") continue;
+       recTime.setTime(timeDirs[timeI], timeI);
+       t = recTime.value();
+       if(t < startTime) continue;
+       if(t > endTime) continue;
+       shift++;
+    }
+
+    scalar dt = origTimeRange / (shift - 1.0);
+    recTime.setEndTime(startTime + 2 * origTimeRange + dt);
+
+    label cellI_transformed = -1;
+    forAll(timeDirs, timeI)
+    {
+       recTime.setTime(timeDirs[timeI], timeI);
+       t = recTime.value();
+       if(t < startTime) continue;
+       if(t > endTime) continue;
+       Info << "time = " << t << ", time index = " << timeI << endl;
+
+       #include "readFields.H"
+
+       forAll(transformedField, cellI)
+       {
+           vector position = mesh.C()[cellI];
+           vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
+           cellI_transformed = mesh.findCell(transformedPosition);
+           if(cellI_transformed < 0)
+           {
+               Info << "Couldn't find transformed cell. Stopping." << endl;
+               return 0;
+           }
+
+           vector value = origField[cellI_transformed];
+           vector transformedValue = -2 * (value & refDirection) * refDirection / (refDirection & refDirection) + value;
+
+           transformedField[cellI] = transformedValue;
+       }
+
+       shiftedTimeI = timeI + shift;
+       t = recTime.value() + origTimeRange + dt;
+       runTime.setTime(t, shiftedTimeI);
+       Info << "creating transformed fields for time = " << t << ", time index = " << shiftedTimeI << endl;
+       transformedField.write();
+    }
+
+    Info << "\nEnd" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirrorVec/readFields.H

@@ -0,0 +1,14 @@
+
+volVectorField origField
+    (
+        IOobject
+        (
+            fieldName,
+            recTime.timePath(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh
+    );
+

applications/utilities/rSmoother/Make/files

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

applications/utilities/rSmoother/Make/options

@@ -0,0 +1,19 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -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)\
+    -lrecurrence \
+    -lfiniteVolume \
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+    -lsampling \
+    -lfvOptions \
+$(CFDEM_ADD_LIBS)

applications/utilities/rSmoother/createFields.H

@@ -0,0 +1,28 @@
+    IOdictionary recProperties
+    (
+        IOobject
+        (
+            "recProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+
+    scalar threshold(readScalar(recProperties.lookup("threshold")));
+
+    volVectorField U_smooth
+    (
+        IOobject
+        (
+            "U_smooth",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0,0,0), vector(0.0, 0.0, 0.0))
+    );