diff --git a/applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C b/applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C
index 54f64497f4b533c52a0b75cab6edb76104674ef3..39774a3373f93ce10a9fafe1ef96d30bdc2671a4 100644
--- a/applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C
+++ b/applications/solvers/incompressible/icoDyMFoam/icoDyMFoam.C
@@ -75,6 +75,9 @@ int main(int argc, char *argv[])
 #           include "correctPhi.H"
         }
 
+        // Keep the absolute fluxes for use in ddtPhiCorr
+        surfaceScalarField phiAbs("phiAbs", phi);
+
         // Make the fluxes relative to the mesh motion
         fvc::makeRelative(phi, U);
 
@@ -95,7 +98,11 @@ int main(int argc, char *argv[])
 
                 U = rAU*UEqn.H();
                 phi = (fvc::interpolate(U) & mesh.Sf());
-                    //+ fvc::ddtPhiCorr(rAU, U, phi);
+
+                if (ddtPhiCorr)
+                {
+                    phi += fvc::ddtPhiCorr(rAU, U, phiAbs);
+                }
 
                 adjustPhi(phi, U, p);
 
@@ -116,7 +123,7 @@ int main(int argc, char *argv[])
                     {
                         pEqn.solve(mesh.solver(p.name()));
                     }
-                
+
                     if (nonOrth == nNonOrthCorr)
                     {
                         phi -= pEqn.flux();
diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
index b3046c7fbb9a55d16fbe357b67995a068b8f2364..e46cba7db9a72752fef3c213accdba95658d49ef 100644
--- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // Make the fluxes absolute
-        if (mesh.changing() && correctPhi)
+        if (mesh.changing())
         {
             fvc::makeAbsolute(phi, U);
         }
@@ -99,8 +99,11 @@ int main(int argc, char *argv[])
             #include "correctPhi.H"
         }
 
+        // Keep the absolute fluxes for use in ddtPhiCorr
+        surfaceScalarField phiAbs("phiAbs", phi);
+
         // Make the fluxes relative to the mesh motion
-        if (mesh.changing() && correctPhi)
+        if (mesh.changing())
         {
             fvc::makeRelative(phi, U);
         }
diff --git a/applications/solvers/multiphase/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interDyMFoam/pEqn.H
index 6ae105b3bcf20acc3a1ece22c980925971bb9098..d4afbae6f49f9f9abd227b2afbb87a7615c4c889 100644
--- a/applications/solvers/multiphase/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interDyMFoam/pEqn.H
@@ -5,12 +5,12 @@
     volVectorField HU = UEqn.H();
     U = rAU*HU;
 
-    surfaceScalarField phiU
-    (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
-    //+ fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
-    );
+    surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));
+
+    if (ddtPhiCorr)
+    {
+        phiU += fvc::ddtPhiCorr(rAU, rho, U, phiAbs);
+    }
 
     phi = phiU +
         (
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
index 61e33a5d7707c52b9d4fc6030330bf741438d316..126b5b428d804c53f4e1d778f14774c8d71615c1 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
@@ -39,6 +39,7 @@ License
 #include "vtkDataArraySelection.h"
 #include "vtkDataSet.h"
 #include "vtkFieldData.h"
+#include "vtkInformation.h"
 #include "vtkMultiBlockDataSet.h"
 #include "vtkRenderer.h"
 #include "vtkTextActor.h"
@@ -61,7 +62,8 @@ void Foam::vtkPV3Foam::AddToBlock
     vtkMultiBlockDataSet* output,
     unsigned int blockNo,
     unsigned int datasetNo,
-    vtkDataSet* dataset
+    vtkDataSet* dataset,
+    const char* datasetName
 )
 {
     vtkDataObject* blockDO = output->GetBlock(blockNo);
@@ -81,6 +83,10 @@ void Foam::vtkPV3Foam::AddToBlock
     }
 
     block->SetBlock(datasetNo, dataset);
+    block->GetMetaData(datasetNo)->Set
+    (
+        vtkCompositeDataSet::NAME(), datasetName
+    );
 }
 
 
@@ -138,25 +144,22 @@ void Foam::vtkPV3Foam::resetCounters()
 }
 
 
-void Foam::vtkPV3Foam::SetName
+void Foam::vtkPV3Foam::SetBlockName
 (
-    vtkUnstructuredGrid* vtkMesh,
+    vtkMultiBlockDataSet* blocks,
+    const int id,
     const char* name
 )
 {
     if (debug)
     {
-        Info<< "entered Foam::vtkPV3Foam::setName" << endl;
-    }
-    vtkCharArray* nmArray =  vtkCharArray::New();
-    nmArray->SetName("Name");
-    size_t len = strlen(name);
-    nmArray->SetNumberOfTuples(static_cast<vtkIdType>(len) + 1);
-    char* copy = nmArray->GetPointer(0);
-    memcpy(copy, name, len);
-    copy[len] = '\0';
-    vtkMesh->GetFieldData()->AddArray(nmArray);
-    nmArray->Delete();
+        Info<< "entered Foam::vtkPV3Foam::setBlockName" << endl;
+    }
+
+    if (blocks->GetMetaData(id) != NULL)
+    {
+        blocks->GetMetaData(id)->Set(vtkCompositeDataSet::NAME(), name);
+    }
 }
 
 
@@ -423,7 +426,7 @@ Foam::vtkPV3Foam::vtkPV3Foam
 
     // Set initial cloud name
     // TODO - TEMPORARY MEASURE UNTIL CAN PROCESS MULTIPLE CLOUDS
-    cloudName_ = "";
+    cloudName_ = "cloud1";
 }
 
 
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
index 8657a9b3e29bffcd88cb1856c46b05cf2bab6937..a2c392bfdad59ead69e017c60e396ce9d41ff515 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
@@ -196,7 +196,8 @@ class vtkPV3Foam
             vtkMultiBlockDataSet* output,
             unsigned int blockNo,
             unsigned int datasetNo,
-            vtkDataSet* dataset
+            vtkDataSet* dataset,
+            const char* datasetName
         );
 
         // Convenience method use to convert the readers from VTK 5
@@ -443,8 +444,13 @@ class vtkPV3Foam
             );
 
 
-        //- Set the name VTK objects
-        void SetName(vtkUnstructuredGrid* vtkMesh, const char* name);
+        //- Set the name of the block
+        void SetBlockName
+        (
+            vtkMultiBlockDataSet* blocks,
+            const int id,
+            const char* name
+        );
 
         //- Disallow default bitwise copy construct
         vtkPV3Foam(const vtkPV3Foam&);
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddFaceSetMesh.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddFaceSetMesh.C
index 9440c2af296612c1bc47dbdb6c004ff19dfa82af..efa6717ec53f0c89bb18bac2ae0174724a6de7db 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddFaceSetMesh.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddFaceSetMesh.C
@@ -51,8 +51,6 @@ void Foam::vtkPV3Foam::addFaceSetMesh
         Info<< "entered add face set internal mesh" << endl;
     }
 
-    SetName(vtkMesh, "faceSetMesh");
-
     // Construct primitivePatch of faces in fSet.
 
     const faceList& faces = mesh.faces();
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddLagrangianMesh.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddLagrangianMesh.C
index c299bd33a90ac937a72fcf6e98691b0d27a38461..dd28a8f930422ee1d2489d30638fcd15cb820f49 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddLagrangianMesh.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddLagrangianMesh.C
@@ -51,8 +51,6 @@ void Foam::vtkPV3Foam::addLagrangianMesh
         Info<< "entered add Lagrangian mesh" << endl;
     }
 
-    SetName(vtkMesh, "LagrangianMesh");
-
     fileNameList cloudDirs
     (
         readDir(mesh.time().timePath()/"lagrangian", fileName::DIRECTORY)
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddPatchMesh.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddPatchMesh.C
index 14e2110a9fe31f88c258c187226c49aa5a793896..08ffc11479ca90dd343bf97c038a3bb16568af34 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddPatchMesh.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddPatchMesh.C
@@ -49,8 +49,6 @@ void Foam::vtkPV3Foam::addPatchMesh
         Info<< "Adding patch: " << p.name() << endl;
     }
 
-    SetName(vtkPatch, p.name().c_str());
-
     if (debug)
     {
         Info<< "converting points" << endl;
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddPointSetMesh.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddPointSetMesh.C
index e28714efeac8ecae3001c65a37116d73e39129c5..26003970ad81a0b8bf04de592e56c3bfe9da105c 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddPointSetMesh.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddPointSetMesh.C
@@ -50,8 +50,6 @@ void Foam::vtkPV3Foam::addPointSetMesh
         Info<< "entered add point set mesh" << endl;
     }
 
-    SetName(vtkMesh, "pointSetMesh");
-
     vtkPoints *vtkpoints = vtkPoints::New();
     vtkpoints->Allocate(mesh.nPoints());
 
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddVolumeMesh.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddVolumeMesh.C
index b1b7a417e8b39b7daab01d268c0ccedf7a114bd6..5d43b29e77b4f6b652c7e2f376ed045cd2b3baf9 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddVolumeMesh.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamAddVolumeMesh.C
@@ -57,8 +57,6 @@ void Foam::vtkPV3Foam::addVolumeMesh
         Info<< "entered add volume mesh" << endl;
     }
 
-    SetName(vtkMesh, "internalMesh");
-
     // Number of additional points needed by the decomposition of polyhedra
     label nAddPoints = 0;
 
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamConvertMesh.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamConvertMesh.C
index 2edace96dcf5fcd903fcd4e69f7fe1417a8e16e6..05ab9b1abfe638e32f79369c754876c15f880fa4 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamConvertMesh.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamConvertMesh.C
@@ -63,10 +63,9 @@ void Foam::vtkPV3Foam::convertMeshVolume
         const fvMesh& mesh = *meshPtr_;
 
         vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::New();
-        SetName(ugrid, "internalMesh");
         addVolumeMesh(mesh, ugrid, superCells_);
-        AddToBlock(output, VOLUME, 0, ugrid);
-//        reader_->SetBlock(output->GetNumberOfBlocks(), ugrid);
+        AddToBlock(output, VOLUME, 0, ugrid, "internalMesh");
+        SetBlockName(output, VOLUME, "Volume");
         selectedRegionDatasetIds_[VOLUME] = 0;
         ugrid->Delete();
     }
@@ -97,7 +96,8 @@ void Foam::vtkPV3Foam::convertMeshLagrangian
 
             vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::New();
             addLagrangianMesh(mesh, ugrid);
-            AddToBlock(output, LAGRANGIAN, 0, ugrid);
+            AddToBlock(output, LAGRANGIAN, 0, ugrid, cloudName_.c_str());
+            SetBlockName(output, LAGRANGIAN, "Lagrangian");
             selectedRegionDatasetIds_[LAGRANGIAN] = 0;
             ugrid->Delete();
         }
@@ -143,7 +143,7 @@ void Foam::vtkPV3Foam::convertMeshPatches
                     .findPatchID(regionName);
                 addPatchMesh(patches[patchId], ugrid);
                 const label nextId = GetNumberOfDataSets(output, VOLUME);
-                AddToBlock(output, VOLUME, nextId, ugrid);
+                AddToBlock(output, VOLUME, nextId, ugrid, regionName.c_str());
                 selectedRegionDatasetIds_[i] = nextId;
                 ugrid->Delete();
             }
@@ -188,7 +188,6 @@ void Foam::vtkPV3Foam::convertMeshCellSet
                 subsetter.setLargeCellSubset(cSet);
 
                 vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::New();
-                SetName(ugrid, cSetName.c_str());
                 const label nextId = GetNumberOfDataSets(output, CELLSET);
                 addVolumeMesh
                 (
@@ -196,11 +195,13 @@ void Foam::vtkPV3Foam::convertMeshCellSet
                     ugrid,
                     superCellSetCells_[nextId]
                 );
-                AddToBlock(output, CELLSET, nextId, ugrid);
+                AddToBlock(output, CELLSET, nextId, ugrid, cSetName.c_str());
                 selectedRegionDatasetIds_[i] = nextId;
                 ugrid->Delete();
             }
         }
+
+        SetBlockName(output, CELLSET, "CellSets");
     }
 }
 
@@ -239,7 +240,6 @@ void Foam::vtkPV3Foam::convertMeshFaceSet
                 const faceSet fSet(mesh, fSetName);
 
                 vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::New();
-                SetName(ugrid, fSetName.c_str());
                 addFaceSetMesh
                 (
                     mesh,
@@ -247,11 +247,13 @@ void Foam::vtkPV3Foam::convertMeshFaceSet
                     ugrid
                 );
                 const label nextId = GetNumberOfDataSets(output, FACESET);
-                AddToBlock(output, FACESET, nextId, ugrid);
+                AddToBlock(output, FACESET, nextId, ugrid, fSetName.c_str());
                 selectedRegionDatasetIds_[i] = nextId;
                 ugrid->Delete();
             }
         }
+
+        SetBlockName(output, FACESET, "FaceSets");
     }
 }
 
@@ -290,7 +292,6 @@ void Foam::vtkPV3Foam::convertMeshPointSet
                 const pointSet pSet(mesh, pSetName);
 
                 vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::New();
-                SetName(ugrid, pSetName.c_str());
                 addPointSetMesh
                 (
                     mesh,
@@ -298,11 +299,13 @@ void Foam::vtkPV3Foam::convertMeshPointSet
                     ugrid
                 );
                 label nextId = GetNumberOfDataSets(output, POINTSET);
-                AddToBlock(output, POINTSET, nextId, ugrid);
+                AddToBlock(output, POINTSET, nextId, ugrid, pSetName.c_str());
                 selectedRegionDatasetIds_[i] = nextId;
                 ugrid->Delete();
             }
         }
+
+        SetBlockName(output, POINTSET, "PointSets");
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/include/readPISOControls.H b/src/finiteVolume/cfdTools/general/include/readPISOControls.H
index d29a60785addc2a2334ef8d1d772aebb45741778..e51ecdbb8bc237b527daf3cbf741b744a4592ab4 100644
--- a/src/finiteVolume/cfdTools/general/include/readPISOControls.H
+++ b/src/finiteVolume/cfdTools/general/include/readPISOControls.H
@@ -25,3 +25,9 @@
     {
         nOuterCorr = readInt(piso.lookup("nOuterCorrectors"));
     }
+
+    bool ddtPhiCorr = false;
+    if (piso.found("ddtPhiCorr"))
+    {
+        ddtPhiCorr = Switch(piso.lookup("ddtPhiCorr"));
+    }
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index 1a5d4b498dd0df79b81914a7ebad728eb15a2eb4..ee3381932e5f0faecf8f0e7f881b681eeb67d737 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+/*---------------------------------------------------------------------------* \
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
@@ -21,7 +21,7 @@ License
     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
-    
+
 \*---------------------------------------------------------------------------*/
 
 #include "EulerDdtScheme.H"
@@ -279,7 +279,7 @@ EulerDdtScheme<Type>::fvmDdt
     scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 
     fvm.diag() = rDeltaT*mesh().V();
-    
+
     if (mesh().moving())
     {
         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
@@ -314,7 +314,7 @@ EulerDdtScheme<Type>::fvmDdt
     scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 
     fvm.diag() = rDeltaT*rho.value()*mesh().V();
-    
+
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
@@ -375,50 +375,30 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
 (
     const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const fluxFieldType& phiAbs
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 
     IOobject ddtIOobject
     (
-        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phiAbs.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rDeltaT.dimensions()*rA.dimensions()*phi.dimensions(),
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        tmp<fluxFieldType> phiCorr =
-            phi.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
+    tmp<fluxFieldType> phiCorr =
+        phiAbs.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
 
-        return tmp<fluxFieldType>
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
         (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr())
-               *fvc::interpolate(rDeltaT*rA)*phiCorr
-            )
-        );
-    }
+            ddtIOobject,
+            fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime(), phiCorr())
+           *fvc::interpolate(rDeltaT*rA)*phiCorr
+        )
+    );
 }
 
 
@@ -429,7 +409,7 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
     const volScalarField& rA,
     const volScalarField& rho,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const fluxFieldType& phiAbs
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
@@ -437,111 +417,94 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
     IOobject ddtIOobject
     (
         "ddtPhiCorr("
-      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+      + rA.name() + ','
+      + rho.name() + ','
+      + U.name() + ','
+      + phiAbs.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    if (mesh().moving())
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phiAbs.dimensions() == dimVelocity*dimArea
+    )
     {
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rDeltaT.dimensions()*rA.dimensions()
-                   *rho.dimensions()*phi.dimensions(),
-                    pTraits<typename flux<Type>::type>::zero
+                rDeltaT
+               *fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
+               *(
+                   fvc::interpolate(rA*rho.oldTime())*phiAbs.oldTime()
+                 - (fvc::interpolate(rA*rho.oldTime()*U.oldTime())
+                  & mesh().Sf())
                 )
             )
         );
     }
-    else
+    else if
+    (
+        U.dimensions() == dimVelocity
+     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
     {
-        if
+        return tmp<fluxFieldType>
         (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
+            new fluxFieldType
             (
-                new fluxFieldType
+                ddtIOobject,
+                rDeltaT
+               *fvcDdtPhiCoeff
                 (
-                    ddtIOobject,
-                    rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
-                   *(
-                        fvc::interpolate(rA*rho.oldTime())*phi.oldTime()
-                      - (fvc::interpolate(rA*rho.oldTime()*U.oldTime())
-                      & mesh().Sf())
-                    )
+                    U.oldTime(),
+                    phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
                 )
-            );
-        }
-        else if 
+               *(
+                   fvc::interpolate(rA*rho.oldTime())
+                  *phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
+                 - (
+                       fvc::interpolate
+                       (
+                           rA*rho.oldTime()*U.oldTime()
+                       ) & mesh().Sf()
+                 )
+               )
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        return tmp<fluxFieldType>
         (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
+            new fluxFieldType
             (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT
-                   *fvcDdtPhiCoeff
-                    (
-                        U.oldTime(),
-                        phi.oldTime()/fvc::interpolate(rho.oldTime())
-                    )
-                   *(
-                        fvc::interpolate(rA*rho.oldTime())
-                       *phi.oldTime()/fvc::interpolate(rho.oldTime())
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*rho.oldTime()*U.oldTime()
-                            ) & mesh().Sf()
-                        )
-                    )
+                ddtIOobject,
+                rDeltaT
+               *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phiAbs.oldTime())
+               *(
+                   fvc::interpolate(rA)*phiAbs.oldTime()
+                 - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
                 )
-            );
-        }
-        else if 
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
         (
-            U.dimensions() == rho.dimensions()*dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT
-                   *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
-                   *(
-                        fvc::interpolate(rA)*phi.oldTime()
-                      - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
-                    )
-                )
-            );
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "EulerDdtScheme<Type>::fvcDdtPhiCorr"
-            )   << "dimensions of phi are not correct"
-                << abort(FatalError);
+            "EulerDdtScheme<Type>::fvcDdtPhiCorr"
+        )   << "dimensions of phiAbs are not correct"
+            << abort(FatalError);
 
-            return fluxFieldType::null();
-        }
+        return fluxFieldType::null();
     }
 }
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index 0c30c641c3fe94aca521c0a8e35ebbb6e88dcc88..f36312b64f30f8591b3a54dc2c3dedc4b72edc4d 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
@@ -21,7 +21,7 @@ License
     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
-    
+
 \*---------------------------------------------------------------------------*/
 
 #include "backwardDdtScheme.H"
@@ -531,58 +531,38 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
         mesh()
     );
 
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rDeltaT.dimensions()*rA.dimensions()*phi.dimensions(),
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        scalar deltaT = deltaT_();
-        scalar deltaT0 = deltaT0_(U);
+    scalar deltaT = deltaT_();
+    scalar deltaT0 = deltaT0_(U);
 
-        scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-        scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-        scalar coefft0  = coefft + coefft00;
+    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    scalar coefft0  = coefft + coefft00;
 
-        return tmp<fluxFieldType>
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
         (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+            ddtIOobject,
+            rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+           *(
+                fvc::interpolate(rA)
                *(
-                    fvc::interpolate(rA)
-                   *(
-                       coefft0*phi.oldTime()
-                     - coefft00*phi.oldTime().oldTime()
-                    )
-                  - (
-                        fvc::interpolate
+                   coefft0*phi.oldTime()
+                 - coefft00*phi.oldTime().oldTime()
+                )
+              - (
+                    fvc::interpolate
+                    (
+                        rA*
                         (
-                            rA*
-                            (
-                                coefft0*U.oldTime()
-                              - coefft00*U.oldTime().oldTime()
-                            )
-                        ) & mesh().Sf()
-                    )
+                            coefft0*U.oldTime()
+                          - coefft00*U.oldTime().oldTime()
+                        )
+                    ) & mesh().Sf()
                 )
             )
-        );
-    }
+        )
+    );
 }
 
 
@@ -593,7 +573,7 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
     const volScalarField& rA,
     const volScalarField& rho,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const fluxFieldType& phiAbs
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
@@ -601,152 +581,134 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
     IOobject ddtIOobject
     (
         "ddtPhiCorr("
-      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+      + rA.name() + ','
+      + rho.name() + ','
+      + U.name() + ','
+      + phiAbs.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    if (mesh().moving())
+    scalar deltaT = deltaT_();
+    scalar deltaT0 = deltaT0_(U);
+
+    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    scalar coefft0  = coefft + coefft00;
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phiAbs.dimensions() == dimVelocity*dimArea
+    )
     {
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rDeltaT.dimensions()*rA.dimensions()
-                   *rho.dimensions()*phi.dimensions(),
-                    pTraits<typename flux<Type>::type>::zero
+                rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
+               *(
+                    coefft0*fvc::interpolate(rA*rho.oldTime())
+                   *phiAbs.oldTime()
+                  - coefft00*fvc::interpolate(rA*rho.oldTime().oldTime())
+                   *phiAbs.oldTime().oldTime()
+                  - (
+                        fvc::interpolate
+                        (
+                            rA*
+                            (
+                                coefft0*rho.oldTime()*U.oldTime()
+                              - coefft00*rho.oldTime().oldTime()
+                               *U.oldTime().oldTime()
+                            )
+                        ) & mesh().Sf()
+                    )
                 )
             )
         );
     }
-    else
+    else if
+    (
+        U.dimensions() == dimVelocity
+     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
     {
-        scalar deltaT = deltaT_();
-        scalar deltaT0 = deltaT0_(U);
-
-        scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-        scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-        scalar coefft0  = coefft + coefft00;
-
-        if
+        return tmp<fluxFieldType>
         (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
+            new fluxFieldType
             (
-                new fluxFieldType
+                ddtIOobject,
+                rDeltaT
+               *fvcDdtPhiCoeff
                 (
-                    ddtIOobject,
-                    rDeltaT*fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
-                   *(
-                        coefft0*fvc::interpolate(rA*rho.oldTime())
-                       *phi.oldTime()
-                      - coefft00*fvc::interpolate(rA*rho.oldTime().oldTime())
-                       *phi.oldTime().oldTime()
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*
-                                (
-                                    coefft0*rho.oldTime()*U.oldTime()
-                                  - coefft00*rho.oldTime().oldTime()
-                                   *U.oldTime().oldTime()
-                                )
-                            ) & mesh().Sf()
-                        )
-                    )
+                    U.oldTime(),
+                    phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
                 )
-            );
-        }
-        else if 
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT
-                   *fvcDdtPhiCoeff
-                    (
-                        U.oldTime(),
-                        phi.oldTime()/fvc::interpolate(rho.oldTime())
-                    )
+               *(
+                    fvc::interpolate(rA*rho.oldTime())
                    *(
-                        fvc::interpolate(rA*rho.oldTime())
-                       *(
-                           coefft0*phi.oldTime()
-                          /fvc::interpolate(rho.oldTime())
-                         - coefft00*phi.oldTime().oldTime()
-                          /fvc::interpolate(rho.oldTime().oldTime())
-                        )
-                      - (
-                            fvc::interpolate
+                       coefft0*phiAbs.oldTime()
+                      /fvc::interpolate(rho.oldTime())
+                     - coefft00*phiAbs.oldTime().oldTime()
+                      /fvc::interpolate(rho.oldTime().oldTime())
+                    )
+                  - (
+                        fvc::interpolate
+                        (
+                            rA*rho.oldTime()*
                             (
-                                rA*rho.oldTime()*
-                                (
-                                    coefft0*U.oldTime()
-                                  - coefft00*U.oldTime().oldTime()
-                                )
-                            ) & mesh().Sf()
-                        )
+                                coefft0*U.oldTime()
+                              - coefft00*U.oldTime().oldTime()
+                            )
+                        ) & mesh().Sf()
                     )
                 )
-            );
-        }
-        else if 
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        return tmp<fluxFieldType>
         (
-            U.dimensions() == rho.dimensions()*dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
+            new fluxFieldType
             (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    rDeltaT
-                   *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
+                ddtIOobject,
+                rDeltaT
+               *fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phiAbs.oldTime())
+               *(
+                    fvc::interpolate(rA)
                    *(
-                        fvc::interpolate(rA)
-                       *(
-                           coefft0*phi.oldTime()
-                         - coefft00*phi.oldTime().oldTime()
-                        )
-                      - (
-                            fvc::interpolate
+                       coefft0*phiAbs.oldTime()
+                     - coefft00*phiAbs.oldTime().oldTime()
+                    )
+                  - (
+                        fvc::interpolate
+                        (
+                            rA*
                             (
-                                rA*
-                                (
-                                    coefft0*U.oldTime()
-                                  - coefft00*U.oldTime().oldTime()
-                                )
-                            ) & mesh().Sf()
-                        )
+                                coefft0*U.oldTime()
+                              - coefft00*U.oldTime().oldTime()
+                            )
+                        ) & mesh().Sf()
                     )
                 )
-            );
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "backwardDdtScheme<Type>::fvcDdtPhiCorr"
-            )   << "dimensions of phi are not correct"
-                << abort(FatalError);
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "backwardDdtScheme<Type>::fvcDdtPhiCorr"
+        )   << "dimensions of phiAbs are not correct"
+            << abort(FatalError);
 
-            return fluxFieldType::null();
-        }
+        return fluxFieldType::null();
     }
 }
 
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
index 63a5f08f916a6749c0a83ebace48cc2aa775795e..0b3efba158ceb8c69558e7649546db2e7910881d 100644
--- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
+++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
@@ -49,7 +49,7 @@ void Foam::cellPointWeight::findTetrahedron
     }
 
     // Initialise closest triangle variables
-    scalar minUVWClose = GREAT;
+    scalar minUVWClose = VGREAT;
     label pointIClose = 0;
     label faceClose = 0;
 
@@ -198,7 +198,7 @@ void Foam::cellPointWeight::findTriangle
     }
 
     // Initialise closest triangle variables
-    scalar minUVClose = GREAT;
+    scalar minUVClose = VGREAT;
     label pointIClose = 0;
 
     // Decompose each face into triangles, making a tet when
diff --git a/tutorials/engineFoam/Allrun b/tutorials/engineFoam/Allrun
index ccab60c0ae5df211bdfbcb6d600d843066a2097b..1e99e22cd4daa34b894c335e3adf20c62e2b4225 100755
--- a/tutorials/engineFoam/Allrun
+++ b/tutorials/engineFoam/Allrun
@@ -13,7 +13,7 @@ runKivaToFoam ()
         echo "kivaToFoam already run on $1: remove log file to run"
     else
         echo "kivaToFoam: converting kiva file"
-        kivaToFoam . $1 -file $2 > $1/log.kivaToFoam 2>&1
+        kivaToFoam -case $1 -file $2 > $1/log.kivaToFoam 2>&1
     fi
 }
 
@@ -23,7 +23,7 @@ restartApplication ()
         echo "$1 already run on $2: remove log file to run"
     else
         echo "Running $1 on $2"
-        $1 . $2 > $2/log-2.$1 2>&1
+        $1 -case $2 > $2/log-2.$1 2>&1
     fi
 }
 
diff --git a/tutorials/engineFoam/kivaTest/constant/combustionProperties b/tutorials/engineFoam/kivaTest/constant/combustionProperties
index ed1bce2ff3030a99fcb310b25a908eebae22cbb6..d6603f9489ac1f2903a31d33fea80cd68fac9d8e 100644
--- a/tutorials/engineFoam/kivaTest/constant/combustionProperties
+++ b/tutorials/engineFoam/kivaTest/constant/combustionProperties
@@ -46,38 +46,38 @@ GuldersCoeffs
 {
     Methane
     {
-        W               W [0 0 0 0 0 0 0] 0.422;
-        eta             eta [0 0 0 0 0 0 0] 0.15;
-        xi              xi [0 0 0 0 0 0 0] 5.18;
-        alpha           alpha [0 0 0 0 0 0 0] 2;
-        beta            beta [0 0 0 0 0 0 0] -0.5;
-        f               f [0 0 0 0 0 0 0] 2.3;
+        W               0.422;
+        eta             0.15;
+        xi              5.18;
+        alpha           2;
+        beta            -0.5;
+        f               2.3;
     }
     Propane
     {
-        W               W [0 0 0 0 0 0 0] 0.446;
-        eta             eta [0 0 0 0 0 0 0] 0.12;
-        xi              xi [0 0 0 0 0 0 0] 4.95;
-        alpha           alpha [0 0 0 0 0 0 0] 1.77;
-        beta            beta [0 0 0 0 0 0 0] -0.2;
-        f               f [0 0 0 0 0 0 0] 2.3;
+        W               0.446;
+        eta             0.12;
+        xi              4.95;
+        alpha           1.77;
+        beta            -0.2;
+        f               2.3;
     }
     IsoOctane
     {
-        W               W [0 0 0 0 0 0 0] 0.4658;
-        eta             eta [0 0 0 0 0 0 0] -0.326;
-        xi              xi [0 0 0 0 0 0 0] 4.48;
-        alpha           alpha [0 0 0 0 0 0 0] 1.56;
-        beta            beta [0 0 0 0 0 0 0] -0.22;
-        f               f [0 0 0 0 0 0 0] 2.3;
+        W               0.4658;
+        eta             -0.326;
+        xi              4.48;
+        alpha           1.56;
+        beta            -0.22;
+        f               2.3;
     }
 }
 
 ignite          yes;
 
-ignitionSites   
+ignitionSites
 (
-    
+
     {
         location        (0.03 0 0.091);
         diameter        0.002;
diff --git a/tutorials/rhoTurbTwinParcelFoam/Allclean b/tutorials/rhoTurbTwinParcelFoam/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..ad62b421fbaa6db3ae7f2e82a016facdb4fe5cfe
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/Allclean
@@ -0,0 +1,16 @@
+#!/bin/sh
+
+currDir=`pwd`
+application=`basename $currDir`
+cases="simplifiedSiwek"
+
+tutorialPath=`dirname $0`/..
+. $tutorialPath/CleanFunctions
+
+wclean $application
+
+for case in $cases
+do
+    cleanCase $case
+done
+
diff --git a/tutorials/rhoTurbTwinParcelFoam/Allrun b/tutorials/rhoTurbTwinParcelFoam/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..f3255d83ad851436e796d7931c0158297c69b5c9
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/Allrun
@@ -0,0 +1,16 @@
+#!/bin/sh
+
+currDir=`pwd`
+application=`basename $currDir`
+cases="simplifiedSiwek"
+
+tutorialPath=`dirname $0`/..
+. $tutorialPath/RunFunctions
+
+compileApplication $currDir $application
+
+for case in $cases
+do
+    runApplication blockMesh $case
+    runApplication $application $case
+done
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/files b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..b11cae03a931a88a5d305111b38cce5ce25e7bed
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/files
@@ -0,0 +1,3 @@
+rhoTurbTwinParcelFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/rhoTurbTwinParcelFoam
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/options b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..4af7133ad332cad337a01c0dad6f21880272033d
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/Make/options
@@ -0,0 +1,21 @@
+EXE_INC = \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels
+
+EXE_LIBS = \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lthermophysicalFunctions \
+    -lbasicThermophysicalModels \
+    -lcombustionThermophysicalModels \
+    -lspecie \
+    -lradiation \
+    -lcompressibleTurbulenceModels
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/UEqn.H b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..4d58a14da89574c12d4162e343f98eba5cf7085f
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/UEqn.H
@@ -0,0 +1,17 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        thermoCloud1.SU1()
+      + kinematicCloud1.SU1()
+      + rho.dimensionedInternalField()*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/createFields.H b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..1191d94b024057dcefd0930ebf96451bf470140c
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/createFields.H
@@ -0,0 +1,84 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicThermo> thermo
+    (
+        basicThermo::New(mesh)
+    );
+
+    volScalarField& p = thermo->p();
+    volScalarField& h = thermo->h();
+    const volScalarField& psi = thermo->psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo->rho()
+    );
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+#   include "compressibleCreatePhi.H"
+
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo()
+        )
+    );
+
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    pointMesh pMesh(mesh);
+    volPointInterpolation vpi(mesh, pMesh);
+
+    Info<< "Constructing thermoCloud1" << endl;
+    basicThermoCloud thermoCloud1
+    (
+        "thermoCloud1",
+        vpi,
+        rho,
+        U,
+        g,
+        thermo()
+    );
+
+    Info<< "Constructing kinematicCloud1" << endl;
+    basicKinematicCloud kinematicCloud1
+    (
+        "kinematicCloud1",
+        vpi,
+        rho,
+        U,
+        thermo().mu(),
+        g
+    );
+
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/hEqn.H b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/hEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..5359c9c2cecae01b52923ff14943cc4a1373d598
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/hEqn.H
@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, h)
+      + fvm::div(phi, h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        DpDt
+      + thermoCloud1.Sh1()
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo->correct();
+}
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/pEqn.H b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..b506245034010d76f0ad0fb87dd22c8b559f5597
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/pEqn.H
@@ -0,0 +1,68 @@
+rho = thermo->rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(thermo->psi())
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam.C b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..0b9edd0b3c8c0349856373914c4d1dc3b7fff851
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam/rhoTurbTwinParcelFoam.C
@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+    rhoTurbFoam
+
+Description
+    Transient solver for compressible, turbulent flow with two thermo-clouds.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicThermo.H"
+#include "compressible/turbulenceModel/turbulenceModel.H"
+
+#include "basicThermoCloud.H"
+#include "basicKinematicCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+
+#   include "createTime.H"
+#   include "createMesh.H"
+#   include "readEnvironmentalProperties.H"
+#   include "createFields.H"
+#   include "readPISOControls.H"
+#   include "initContinuityErrs.H"
+#   include "readTimeControls.H"
+#   include "compressibleCourantNo.H"
+#   include "setInitialDeltaT.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+#       include "readTimeControls.H"
+#       include "readPISOControls.H"
+#       include "compressibleCourantNo.H"
+#       include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        Info<< "Evolving thermoCloud1" << endl;
+        thermoCloud1.evolve();
+        thermoCloud1.info();
+
+        Info<< "Evolving kinematicCloud1" << endl;
+        kinematicCloud1.evolve();
+        kinematicCloud1.info();
+
+
+#       include "rhoEqn.H"
+
+        // --- PIMPLE loop
+        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+        {
+#           include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=1; corr<=nCorr; corr++)
+            {
+#               include "hEqn.H"
+#               include "pEqn.H"
+            }
+        }
+
+        turbulence->correct();
+
+        rho = thermo->rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/G b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/G
new file mode 100644
index 0000000000000000000000000000000000000000..473b2202225297ae58694e6f19ceced3bdbd8800
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/G
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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           volScalarField;
+    object          G;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+dimensions      [1 0 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    top
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1.0;
+        value           uniform 0;
+    }
+    bottom
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1.0;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1.0;
+        value           uniform 0;
+    }
+    symmetry
+    {
+        type            symmetryPlane;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/T b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/T
new file mode 100644
index 0000000000000000000000000000000000000000..48baa8f7f74f2fc15c44e02ec39a9d49ec075cd7
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/T
@@ -0,0 +1,61 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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           volScalarField;
+    object          T;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 400;
+
+boundaryField
+{
+    top
+    {
+        type            fixedValue;
+        value           uniform 400;
+    }
+
+    bottom
+    {
+        type            zeroGradient;
+    }
+
+    walls
+    {
+        type            zeroGradient;
+    }
+
+    symmetry
+    {
+        type            symmetryPlane;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/U b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..7d433912331103f9cd7544667b4392e5361d9058
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/U
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4.1                                 |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version 2.0;
+    format ascii;
+
+    root "";
+    case "";
+    instance "";
+    local "";
+
+    class volVectorField;
+    object U;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    top
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    bottom
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    walls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    symmetry
+    {
+        type            symmetryPlane;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/epsilon b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..8de558dbd3a81171122c5b74b51a6efaf0604abd
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/epsilon
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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           volScalarField;
+    object          epsilon;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 5390.5;
+
+boundaryField
+{
+    top
+    {
+        type            zeroGradient;
+    }
+
+    bottom
+    {
+        type            zeroGradient;
+    }
+
+    walls
+    {
+        type            zeroGradient;
+    }
+
+    symmetry
+    {
+        type            symmetryPlane;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/k b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..f5f4cb92dd5a03af12d8f68cfe07762ff3b7e90e
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/k
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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           volScalarField;
+    object          k;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 37.5;
+
+boundaryField
+{
+    top
+    {
+        type            zeroGradient;
+    }
+
+    bottom
+    {
+        type            zeroGradient;
+    }
+
+    walls
+    {
+        type            zeroGradient;
+    }
+
+    symmetry
+    {
+        type            symmetryPlane;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/p b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..e2f41ffd04e6bc2aa692c3e4e6b0acd304dc9739
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/0/p
@@ -0,0 +1,2558 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4.1                                 |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version 2.0;
+    format ascii;
+
+    root "/home/andy/OpenFOAM/andy-1.4.1/development/spray/rhoTurbThermoParcelExplicitSourceFoam";
+    case "testCase";
+    instance "0";
+    local "";
+
+    class volScalarField;
+    object p;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   nonuniform List<scalar>
+2500
+(
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+5e+05
+)
+;
+
+boundaryField
+{
+    top
+    {
+        type            zeroGradient;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    symmetry
+    {
+        type            symmetryPlane;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/environmentalProperties b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/environmentalProperties
new file mode 100644
index 0000000000000000000000000000000000000000..e9aee6c9b50caa4fae9f0025596b7894773d2cbb
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/environmentalProperties
@@ -0,0 +1,28 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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          environmentalProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+g               g [0 1 -2 0 0 0 0] (0 -9.81 0);
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Positions b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Positions
new file mode 100644
index 0000000000000000000000000000000000000000..f01a5a8ad8f65ff5c64b2247c279212bf9f9ca2b
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Positions
@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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           vectorField;
+    object          kinematicCloud1Positions;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+(
+(0.0075 0.5 0.05)
+(0.0125 0.5 0.05)
+(0.0175 0.5 0.05)
+(0.0225 0.5 0.05)
+(0.0275 0.5 0.05)
+(0.0325 0.5 0.05)
+(0.0375 0.5 0.05)
+(0.0425 0.5 0.05)
+(0.0475 0.5 0.05)
+(0.0075 0.4 0.05)
+(0.0125 0.4 0.05)
+(0.0175 0.4 0.05)
+(0.0225 0.4 0.05)
+(0.0275 0.4 0.05)
+(0.0325 0.4 0.05)
+(0.0375 0.4 0.05)
+(0.0425 0.4 0.05)
+(0.0475 0.4 0.05)
+)
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
new file mode 100644
index 0000000000000000000000000000000000000000..3d652853607546d2299e6b1cfbd0eedbb8094897
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4.1                                 |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          kinematicCloud1Properties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Particle sub-models
+InjectionModel                           ManualInjection;
+DragModel                                SphereDrag;
+DispersionModel                          StochasticDispersionRAS;
+WallInteractionModel                     StandardWallInteraction;
+
+// Parcel basis type
+parcelBasisType                          mass;
+
+// Total mass to inject
+massTotal  massTotal [ 1  0  0  0  0]    2.0e-4;
+
+// Minimum particle mass
+minParticleMass      minParticleMass     [ 1  0  0  0  0]     1.0e-15;
+
+// Parcel thermo properties
+rho0      rho0     [ 1 -3  0  0  0]      5000;
+
+// Coupling between particles and carrier phase via source terms
+coupled                                  true;
+
+// Integer used to identify different parcel types
+parcelTypeId                             2;
+
+interpolationSchemes
+{
+    rho                                  cell;
+    U                                    cellPointFace;
+    mu                                   cell;
+}
+
+integrationSchemes
+{
+    U                                    Euler;
+}
+
+ManualInjectionCoeffs
+{
+    injectionTime                        0;
+    positionsFile                        kinematicCloud1Positions;
+    U0                                   (0 0 0);
+    parcelPDF
+    {
+        pdfType                          RosinRammler;
+        RosinRammlerPDF
+        {
+            minValue                     50.0e-06;
+            maxValue                     100.0e-06;
+            d                            (75.0e-06);
+            n                            (0.5);
+        }
+    }
+}
+
+StandardWallInteractionCoeffs
+{
+    e      e        [ 0  0  0  0  0]     1;
+    mu     mu       [ 0  0  0  0  0]     0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/polyMesh/blockMeshDict b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..26ea6f0982bf65f1dd8c723e8dfaa507bff37fdf
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/polyMesh/blockMeshDict
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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          blockMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1.0;
+
+vertices
+(
+    (0     0     0)
+    (0.05  0     0)
+    (0.05  0.5   0)
+    (0     0.5   0)
+    (0     0     0.1)
+    (0.05  0     0.1)
+    (0.05  0.5   0.1)
+    (0     0.5   0.1)
+    (0.5   0     0)
+    (0.5   0.5   0)
+    (0.5   0     0.1)
+    (0.5   0.5   0.1)
+    (0.05  1     0)
+    (0     1     0)
+    (0.05  1     0.1)
+    (0     1     0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7)     (5 50 1) simpleGrading (1 1 1)
+    hex (1 8 9 2 5 10 11 6)   (40 50 1) simpleGrading (1 1 1)
+    hex (3 2 12 13 7 6 14 15) (5 50 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    patch top
+    (
+        (13 15 14 12)
+    )
+    patch bottom
+    (
+        (0 1 5 4)
+        (1 8 10 5)
+    )
+    wall walls
+    (
+        (8 9 11 10)
+        (9 2 6 11)
+        (2 12 14 6)
+    )
+    symmetryPlane symmetry
+    (
+        (4 7 3 0)
+        (7 15 13 3)
+    )
+    empty frontAndBack
+    (
+        (0 3 2 1)
+        (3 13 12 2)
+        (1 2 9 8)
+        (5 6 7 4)
+        (6 14 15 7)
+        (10 11 6 5)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Positions b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Positions
new file mode 100644
index 0000000000000000000000000000000000000000..ced7be69f3489984f0ce419225dc61a4bd55c5f1
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Positions
@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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           vectorField;
+    object          limestonePositions;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+(
+(0.0075 0.55 0.05)
+(0.0125 0.55 0.05)
+(0.0175 0.55 0.05)
+(0.0225 0.55 0.05)
+(0.0275 0.55 0.05)
+(0.0325 0.55 0.05)
+(0.0375 0.55 0.05)
+(0.0425 0.55 0.05)
+(0.0475 0.55 0.05)
+(0.0075 0.45 0.05)
+(0.0125 0.45 0.05)
+(0.0175 0.45 0.05)
+(0.0225 0.45 0.05)
+(0.0275 0.45 0.05)
+(0.0325 0.45 0.05)
+(0.0375 0.45 0.05)
+(0.0425 0.45 0.05)
+(0.0475 0.45 0.05)
+)
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
new file mode 100644
index 0000000000000000000000000000000000000000..af6f908566aeb3e8bc53d98e05f646c9dc5e1896
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4.1                                 |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          thermoCloud1Properties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Particle sub-models
+InjectionModel                           ManualInjection;
+DragModel                                SphereDrag;
+DispersionModel                          StochasticDispersionRAS;
+WallInteractionModel                     StandardWallInteraction;
+HeatTransferModel                        RanzMarshall;
+
+radiation                                off;
+
+// Parcel basis type
+parcelBasisType                          mass;
+
+// Total mass to inject
+massTotal  massTotal [ 1  0  0  0  0]    1e-4;
+
+// Minimum particle mass
+minParticleMass      minParticleMass     [ 1  0  0  0  0]     1.0e-15;
+
+// Parcel thermo properties
+rho0      rho0     [ 1 -3  0  0  0]      2500;
+T0        T0       [ 0  0  0  1  0]      300;
+cp0       cp0      [ 0  2 -2 -1  0]      900;
+epsilon0  epsilon0 [ 0  0  0  0  0]      1;
+f0        f0       [ 0  0  0  0  0]      0.5;
+
+// Coupling between particles and carrier phase via source terms
+coupled                                  true;
+
+// Integer used to identify different parcel types
+parcelTypeId                             1;
+
+interpolationSchemes
+{
+    rho                                  cell;
+    U                                    cellPointFace;
+    mu                                   cell;
+    T                                    cell;
+    Cp                                   cell;
+}
+
+integrationSchemes
+{
+    U                                    Euler;
+    T                                    Analytical;
+}
+
+ManualInjectionCoeffs
+{
+    injectionTime                        0;
+    positionsFile                        thermoCloud1Positions;
+    U0                                   (0 0 0);
+    parcelPDF
+    {
+        pdfType                          RosinRammler;
+        RosinRammlerPDF
+        {
+            minValue                     5.0e-06;
+            maxValue                     500.0e-06;
+            d                            (50.0e-06);
+            n                            (0.5);
+        }
+    }
+}
+
+StandardWallInteractionCoeffs
+{
+    e      e        [ 0  0  0  0  0]     1;
+    mu     mu       [ 0  0  0  0  0]     0;
+}
+
+RanzMarshallCoeffs
+{
+    Pr    Pr       [ 0  0  0  0  0]     0.7;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermophysicalProperties b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermophysicalProperties
new file mode 100644
index 0000000000000000000000000000000000000000..de5bb8e48e5192cafcbc12ff31c5be7686645f48
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/thermophysicalProperties
@@ -0,0 +1,32 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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          thermophysicalProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Thermophysical model
+thermoType      hThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+
+mixture         air 1 28.9 1007 0 1.84e-05 0.7;
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/turbulenceProperties b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..eb33720fe872b2016d815aea5fd8fd6634f22155
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/constant/turbulenceProperties
@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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          turbulenceProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Turbulence model selection
+turbulenceModel kEpsilon;
+
+// Do you wish to calculate turbulence?
+turbulence      on;
+
+// Laminar model coefficients
+laminarCoeffs
+{
+}
+
+// Standard k-epsilon model coefficients
+kEpsilonCoeffs
+{
+    // Cmu
+    Cmu             Cmu [0 0 0 0 0 0 0] 0.09;
+    // C1
+    C1              C1 [0 0 0 0 0 0 0] 1.44;
+    // C2
+    C2              C2 [0 0 0 0 0 0 0] 1.92;
+    // C3
+    C3              C3 [0 0 0 0 0 0 0] -0.33;
+    // alphah
+    alphah          alphah [0 0 0 0 0 0 0] 1;
+    // alphak
+    alphak          alphak [0 0 0 0 0 0 0] 1;
+    // alphaEps
+    alphaEps        alphaEps [0 0 0 0 0 0 0] 0.76923;
+}
+
+// RNG k-epsilon model coefficients
+RNGkEpsilonCoeffs
+{
+    // Cmu
+    Cmu             Cmu [0 0 0 0 0 0 0] 0.0845;
+    // C1
+    C1              C1 [0 0 0 0 0 0 0] 1.42;
+    // C2
+    C2              C2 [0 0 0 0 0 0 0] 1.68;
+    // C3
+    C3              C3 [0 0 0 0 0 0 0] -0.33;
+    // alphah
+    alphah          alphah [0 0 0 0 0 0 0] 1;
+    // alphak
+    alphak          alphaK [0 0 0 0 0 0 0] 1.39;
+    // alphaEps
+    alphaEps        alphaEps [0 0 0 0 0 0 0] 1.39;
+    // eta0
+    eta0            eta0 [0 0 0 0 0 0 0] 4.38;
+    // beta
+    beta            beta [0 0 0 0 0 0 0] 0.012;
+}
+
+// Launder-Sharma low Reynolds number k-epsilon model coefficients
+LaunderSharmaKECoeffs
+{
+    // Cmu
+    Cmu             Cmu [0 0 0 0 0 0 0] 0.09;
+    // C1
+    C1              C1 [0 0 0 0 0 0 0] 1.44;
+    // C2
+    C2              C2 [0 0 0 0 0 0 0] 1.92;
+    // C3
+    C3              C3 [0 0 0 0 0 0 0] -0.33;
+    // alphah
+    alphah          alphah [0 0 0 0 0 0 0] 1;
+    // alphak
+    alphak          alphak [0 0 0 0 0 0 0] 1;
+    // alphaEps
+    alphaEps        alphaEps [0 0 0 0 0 0 0] 0.76923;
+}
+
+// Launder-Reece-Rodi RSTM with wall functions model coefficients
+LRRCoeffs
+{
+    // Cmu
+    Cmu             Cmu [0 0 0 0 0 0 0] 0.09;
+    // Clrr1
+    Clrr1           Clrr1 [0 0 0 0 0 0 0] 1.8;
+    // Clrr2
+    Clrr2           Clrr2 [0 0 0 0 0 0 0] 0.6;
+    // C1
+    C1              C1 [0 0 0 0 0 0 0] 1.44;
+    // C2
+    C2              C2 [0 0 0 0 0 0 0] 1.92;
+    // Cs
+    Cs              Cs [0 0 0 0 0 0 0] 0.25;
+    // Ceps
+    Ceps            Ceps [0 0 0 0 0 0 0] 0.15;
+    // alphah
+    alphah          alphah [0 0 0 0 0 0 0] 1;
+    // alphaEps
+    alphaEps        alphaEps [0 0 0 0 0 0 0] 0.76923;
+    // alphaR
+    alphaR          alphaR [0 0 0 0 0 0 0] 1.22;
+}
+
+// Launder-Gibson RSTM with wall reflection and wall functions model coefficients
+LaunderGibsonRSTMCoeffs
+{
+    // Cmu
+    Cmu             Cmu [0 0 0 0 0 0 0] 0.09;
+    // Clg1
+    Clg1            Clg1 [0 0 0 0 0 0 0] 1.8;
+    // Clg2
+    Clg2            Clg2 [0 0 0 0 0 0 0] 0.6;
+    // C1
+    C1              C1 [0 0 0 0 0 0 0] 1.44;
+    // C2
+    C2              C2 [0 0 0 0 0 0 0] 1.92;
+    // C1Ref
+    C1Ref           C1Ref [0 0 0 0 0 0 0] 0.5;
+    // C2Ref
+    C2Ref           C2Ref [0 0 0 0 0 0 0] 0.3;
+    // Cs
+    Cs              Cs [0 0 0 0 0 0 0] 0.25;
+    // Ceps
+    Ceps            Ceps [0 0 0 0 0 0 0] 0.15;
+    // alphah
+    alphah          alphah [0 0 0 0 0 0 0] 1;
+    // alphaEps
+    alphaEps        alphaEps [0 0 0 0 0 0 0] 0.76923;
+    // alphaR
+    alphaR          alphaR [0 0 0 0 0 0 0] 1.22;
+}
+
+// Wall function coefficients
+wallFunctionCoeffs
+{
+    // kappa
+    kappa           kappa [0 0 0 0 0 0 0] 0.4187;
+    // E
+    E               E [0 0 0 0 0 0 0] 9;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/controlDict b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..05d2cad84f04614e8e64bafe2305d7230ed42f7f
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/controlDict
@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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          controlDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Foam Application Class
+application     rhoTurbThermoParcelFoam;
+
+// Start point of run
+startFrom       latestTime;
+
+// Calculation start time
+startTime       0;
+
+// End point of run
+stopAt          endTime;
+
+// Calculation end time
+endTime         0.5;
+
+// Calculation time step
+deltaT          1.0e-4;
+
+// Type of write output control
+writeControl    adjustableRunTime;
+
+// Interval with which the results are output
+writeInterval   0.01;
+
+// Limits number of time directories before overwriting
+purgeWrite      0;
+
+// Write Format
+writeFormat     ascii;
+
+// Significant figures of written ASCII data
+writePrecision  10;
+
+// Write Compression
+writeCompression uncompressed;
+
+// Time directories name format
+timeFormat      general;
+
+// Decimal precision of time directory names
+timePrecision   6;
+
+// Can parameters be modified during run time?
+runTimeModifiable yes;
+
+// Automatic adjustment of time step?
+adjustTimeStep  yes;
+
+// maxCo
+maxCo           0.2;
+
+// maxDeltaT
+maxDeltaT       1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/decomposeParDict b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..dea9c89534ea9f14d56b019197dcd1188a669726
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/decomposeParDict
@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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          decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+numberOfSubdomains 4;
+
+method          metis;
+
+simpleCoeffs
+{
+    n               (2 2 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (1 1 1);
+    delta           0.001;
+    order           xyz;
+}
+
+metisCoeffs
+{
+    processorWeights 
+    (
+        1
+        1
+        1
+        1
+    );
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           
+(
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/fvSchemes b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/fvSchemes
new file mode 100755
index 0000000000000000000000000000000000000000..029fee22959778296efa7bc9f2d3261b653198d1
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/fvSchemes
@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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          fvSchemes;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Time derivative discretisation schemes
+ddtSchemes
+{
+    // Default scheme
+    default         Euler;
+}
+
+// Gradient discretisation schemes
+gradSchemes
+{
+    // Default gradient scheme
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+// Convection discretisation schemes
+divSchemes
+{
+    // Default scheme
+    default         none;
+    div(phi,U)      Gauss upwind;
+    div(phid,p)     Gauss upwind;
+    div(phiU,p)     Gauss linear;
+    div(phi,h)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) Gauss upwind;
+    div(U)          Gauss linear;
+    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div(phi,Yi_h)   Gauss  upwind;
+}
+
+// Laplacian discretisation schemes
+laplacianSchemes
+{
+    // Default scheme
+    default                 Gauss linear corrected;
+    laplacian(muEff,U) Gauss linear corrected;
+    laplacian(mut,U) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DREff,R) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(alphaEff,h) Gauss linear corrected;
+}
+
+// Interpolation schemes
+interpolationSchemes
+{
+    // Default scheme
+    default         linear;
+}
+
+// Surface normal gradient schemes
+snGradSchemes
+{
+    // Default scheme
+    default         corrected;
+}
+
+// Calculation of flux
+fluxRequired
+{
+    // Create storage for flux for all solved variables?
+    default         no;
+    p;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/fvSolution b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/fvSolution
new file mode 100755
index 0000000000000000000000000000000000000000..1fc5059b5c248063f54af762a7e9aeff133cc21d
--- /dev/null
+++ b/tutorials/rhoTurbTwinParcelFoam/simplifiedSiwek/system/fvSolution
@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  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          fvSolution;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    // Solver for the rho equation
+    rho PCG
+    {
+        preconditioner   DIC;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    // Solver for the U equation
+    U PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    // Solver for the p equation
+    p PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+    G PCG
+    {
+        preconditioner   DIC;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    Yi PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+    CO2 PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+    O2 PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+    N2 PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+    CH4 PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+    H2 PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+    H2O PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+    CO PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-06;
+        relTol           0;
+    };
+
+    // Solver for the h equation
+    h PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    // Solver for the R equation
+    R PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    // Solver for the k equation
+    k PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    // Solver for the epsilon equation
+    epsilon PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+}
+
+PISO
+{
+    // Transonic?
+    transonic yes;
+    // Number of PISO correctors
+    nCorrectors     2;
+    // Number of non-orthogonal correctors
+    nNonOrthogonalCorrectors 0;
+    // momentumPredictor?
+    momentumPredictor yes;
+}
+
+
+// ************************************************************************* //