WIP: TUT: tutorial adjustment for backward compatibility

WIP: TUT: tutorial adjustments for backward compatibility
WIP: Updated lagrangian solvers followinf library changes
2023-01-31 16:14:43 +00:00 · 2023-01-20 16:56:11 +00:00 · 2023-01-10 13:35:42 +00:00 · 2023-01-10 13:35:42 +00:00 · 2023-01-10 13:35:42 +00:00 · 2023-01-10 13:35:42 +00:00
331 changed files with 4369 additions and 2792 deletions
--- a/META-INFO/api-info
+++ b/META-INFO/api-info
@ -1,2 +1,2 @@
 api=2212
-patch=240625
+patch=0
--- a/applications/solvers/combustion/fireFoam/createClouds.H
+++ b/applications/solvers/combustion/fireFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingCloud parcels
+reactingCloud parcels
 (
    "reactingCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@ -37,7 +37,7 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicReactingCloud.H"
+#include "reactingCloud.H"
 #include "surfaceFilmModel.H"
 #include "pyrolysisModelCollection.H"
 #include "radiationModel.H"
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@ -9,7 +9,7 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@ -7,7 +7,7 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
@ -10,7 +10,7 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
@ -7,7 +7,7 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@ -19,7 +19,7 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude


 EXE_LIBS = \
--- a/applications/solvers/incompressible/pimpleFoam/Make/options
+++ b/applications/solvers/incompressible/pimpleFoam/Make/options
@ -8,7 +8,7 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
@ -42,11 +42,11 @@ Description
 #include "CorrectPhi.H"

 #ifdef MPPIC
-    #include "basicKinematicCloud.H"
-    #define basicKinematicTypeCloud basicKinematicCloud
+    #include "kinematicCloud.H"
+    #define kinematicTypeCloud kinematicCloud
 #else
-    #include "basicKinematicCollidingCloud.H"
-    #define basicKinematicTypeCloud basicKinematicCollidingCloud
+    #include "kinematicCollidingCloud.H"
+    #define kinematicTypeCloud kinematicCollidingCloud
 #endif

 int main(int argc, char *argv[])
@ -88,7 +88,7 @@ int main(int argc, char *argv[])
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // Store the particle positions
-        kinematicCloud.storeGlobalPositions();
+        kCloud.storeGlobalPositions();

        mesh.update();

@ -111,15 +111,15 @@ int main(int argc, char *argv[])
        continuousPhaseTransport.correct();
        muc = rhoc*continuousPhaseTransport.nu();

-        kinematicCloud.evolve();
+        kCloud.evolve();

        // Update continuous phase volume fraction field
-        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+        alphac = max(1.0 - kCloud.theta(), alphacMin);
        alphac.correctBoundaryConditions();
        alphacf = fvc::interpolate(alphac);
        alphaPhic = alphacf*phic;

-        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
+        fvVectorMatrix cloudSU(kCloud.SU(Uc));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
@ -12,6 +12,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
@ -11,8 +11,11 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
--- a/applications/solvers/lagrangian/DPMFoam/DPMFoam.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMFoam.C
@ -43,11 +43,11 @@ Description
 #include "pimpleControl.H"

 #ifdef MPPIC
-    #include "basicKinematicCloud.H"
-    #define basicKinematicTypeCloud basicKinematicCloud
+    #include "kinematicCloud.H"
+    #define kinematicTypeCloud kinematicCloud
 #else
-    #include "basicKinematicCollidingCloud.H"
-    #define basicKinematicTypeCloud basicKinematicCollidingCloud
+    #include "kinematicCollidingCloud.H"
+    #define kinematicTypeCloud kinematicCollidingCloud
 #endif

 int main(int argc, char *argv[])
@ -91,16 +91,16 @@ int main(int argc, char *argv[])
        continuousPhaseTransport.correct();
        muc = rhoc*continuousPhaseTransport.nu();

-        Info<< "Evolving " << kinematicCloud.name() << endl;
-        kinematicCloud.evolve();
+        Info<< "Evolving " << kCloud.name() << endl;
+        kCloud.evolve();

        // Update continuous phase volume fraction field
-        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+        alphac = max(1.0 - kCloud.theta(), alphacMin);
        alphac.correctBoundaryConditions();
        alphacf = fvc::interpolate(alphac);
        alphaPhic = alphacf*phic;

-        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
+        fvVectorMatrix cloudSU(kCloud.SU(Uc));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
--- a/applications/solvers/lagrangian/DPMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/Make/options
@ -10,8 +10,11 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
--- a/applications/solvers/lagrangian/DPMFoam/createFields.H
+++ b/applications/solvers/lagrangian/DPMFoam/createFields.H
@ -125,13 +125,13 @@ const word kinematicCloudName
 );

 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-basicKinematicTypeCloud kinematicCloud
+kinematicTypeCloud kCloud
 (
    kinematicCloudName,
+    g,
    rhoc,
    Uc,
-    muc,
-    g
+    muc
 );

 // Particle fraction upper limit
@ -139,13 +139,13 @@ scalar alphacMin
 (
    1.0
  - (
-        kinematicCloud.particleProperties().subDict("constantProperties")
+        kCloud.particleProperties().subDict("constantProperties")
       .get<scalar>("alphaMax")
    )
 );

 // Update alphac from the particle locations
-alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+alphac = max(1.0 - kCloud.theta(), alphacMin);
 alphac.correctBoundaryConditions();

 surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
--- a/applications/solvers/lagrangian/coalChemistryFoam/Make/options
+++ b/applications/solvers/lagrangian/coalChemistryFoam/Make/options
@ -6,7 +6,6 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@ -35,7 +34,6 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -lcoalCombustion\
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@ -37,8 +37,8 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicThermoCloud.H"
-#include "coalCloud.H"
+#include "thermoCloud.H"
+#include "reactingMultiphaseCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "fvOptions.H"
--- a/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H
@ -1,19 +1,19 @@
 Info<< "\nConstructing coal cloud" << endl;
-coalCloud coalParcels
+reactingMultiphaseCloud coalParcels
 (
    "coalCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );

 Info<< "\nConstructing limestone cloud" << endl;
-basicThermoCloud limestoneParcels
+thermoCloud limestoneParcels
 (
    "limestoneCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/createFields.H
@ -63,13 +63,13 @@ const word kinematicCloudName
 );

 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-basicKinematicCollidingCloud kinematicCloud
+kinematicCollidingCloud kCloud
 (
    kinematicCloudName,
+    g,
    rhoInf,
    U,
-    mu,
-    g
+    mu
 );

 IOobject Hheader
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C
@ -41,7 +41,7 @@ Description
 #include "dynamicFvMesh.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
-#include "basicKinematicCollidingCloud.H"
+#include "kinematicCollidingCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -76,19 +76,19 @@ int main(int argc, char *argv[])
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

-        kinematicCloud.storeGlobalPositions();
+        kCloud.storeGlobalPositions();

        mesh.update();

        U.correctBoundaryConditions();

-        Info<< "Evolving " << kinematicCloud.name() << endl;
+        Info<< "Evolving " << kCloud.name() << endl;

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

-        kinematicCloud.evolve();
+        kCloud.evolve();

        runTime.write();

--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C
@ -40,7 +40,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
-#include "basicKinematicCollidingCloud.H"
+#include "kinematicCollidingCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -75,13 +75,13 @@ int main(int argc, char *argv[])
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

-        Info<< "Evolving " << kinematicCloud.name() << endl;
+        Info<< "Evolving " << kCloud.name() << endl;

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

-        kinematicCloud.evolve();
+        kCloud.evolve();

        runTime.write();

--- a/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
@ -8,6 +8,12 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
--- a/applications/solvers/lagrangian/kinematicParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/UEqn.H
@ -6,7 +6,7 @@
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
-        parcels.SU(U, true)
+        scalar(1)/rhoInfValue*parcels.SU(U)
      + fvOptions(U)
    );

--- a/applications/solvers/lagrangian/kinematicParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/createClouds.H
@ -4,13 +4,12 @@ const word kinematicCloudName
 );

 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-
-basicKinematicCloud parcels
+kinematicCloud parcels
 (
    kinematicCloudName,
+    g,
    rhoInf,
    U,
-    muc,
-    g
+    muc
 );

--- a/applications/solvers/lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
@ -40,7 +40,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "surfaceFilmModel.H"
-#include "basicKinematicCloud.H"
+#include "kinematicCloud.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
 #include "CorrectPhi.H"
@ -67,7 +67,6 @@ int main(int argc, char *argv[])
    #include "createDyMControls.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
-    #include "createRegionControls.H"
    #include "createUfIfPresent.H"

    turbulence->validate();
@ -95,7 +94,7 @@ int main(int argc, char *argv[])
        // Do any mesh changes
        mesh.update();

-        if (solvePrimaryRegion && mesh.changing())
+        if (pimple.solveFlow() && mesh.changing())
        {
            MRF.update();

@ -120,7 +119,7 @@ int main(int argc, char *argv[])
        parcels.evolve();
        surfaceFilm.evolve();

-        if (solvePrimaryRegion)
+        if (pimple.solveFlow())
        {
            // --- PIMPLE loop
            while (pimple.loop())
--- a/applications/solvers/lagrangian/parcelFoam/Make/files
+++ b/applications/solvers/lagrangian/parcelFoam/Make/files
@ -0,0 +1,3 @@
+parcelFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/parcelFoam
--- a/applications/solvers/lagrangian/parcelFoam/Make/options
+++ b/applications/solvers/lagrangian/parcelFoam/Make/options
@ -1,49 +1,48 @@
 EXE_INC = \
+    -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/faOptions/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude

-LIB_LIBS = \
+
+EXE_LIBS = \
    -lfiniteVolume \
+    -lfvOptions \
    -lmeshTools \
+    -lsampling \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -latmosphericModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lsurfaceFilmDerivedFvPatchFields \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -ldistributionModels \
-    -lspecie \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lthermophysicalProperties \
-    -lreactionThermophysicalModels \
-    -lSLGThermo \
-    -lradiationModels \
-    -lturbulenceModels \
-    -lincompressibleTurbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lincompressibleTransportModels \
-    -lregionModels \
-    -lsurfaceFilmModels \
-    -ldynamicMesh \
-    -ldynamicFvMesh \
    -lregionFaModels \
-    -lfiniteArea
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/parcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/parcelFoam/UEqn.H
@ -0,0 +1,22 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(U) + fvm::div(phi, U)
+      + MRF.DDt(U)
+      + turbulence->divDevReff(U)
+     ==
+        scalar(1)/rhoInfValue*parcels.SU(U)
+      + fvOptions(U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve(UEqn == -fvc::grad(p));
+
+        fvOptions.correct(U);
+    }
--- a/applications/solvers/lagrangian/parcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/parcelFoam/createClouds.H
@ -0,0 +1,2 @@
+auto parcels = parcelCloudModelList(g, rhoInf, U, muc);
+
--- a/applications/solvers/lagrangian/parcelFoam/createFieldRefs.H
+++ b/applications/solvers/lagrangian/parcelFoam/createFieldRefs.H
@ -0,0 +1 @@
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
--- a/applications/solvers/lagrangian/parcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/parcelFoam/createFields.H
@ -0,0 +1,85 @@
+#include "readGravitationalAcceleration.H"
+
+Info<< "Reading field p\n" << endl;
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "createPhi.H"
+
+singlePhaseTransportModel laminarTransport(U, phi);
+
+dimensionedScalar rhoInfValue
+(
+    "rhoInf",
+    dimDensity,
+    laminarTransport
+);
+
+volScalarField rhoInf
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    rhoInfValue
+);
+
+volScalarField muc
+(
+    IOobject
+    (
+        "muc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    rhoInf*laminarTransport.nu()
+);
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, laminarTransport)
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+mesh.setFluxRequired(p.name());
+
+#include "createMRF.H"
+#include "createClouds.H"
+#include "createSurfaceFilmModel.H"
+#include "createFvOptions.H"
+
+
--- a/applications/solvers/lagrangian/parcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/parcelFoam/pEqn.H
@ -0,0 +1,58 @@
+
+volScalarField rAU(1.0/UEqn.A());
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
+
+if (pimple.ddtCorr())
+{
+    phiHbyA += MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf));
+}
+
+MRF.makeRelative(phiHbyA);
+
+if (p.needReference())
+{
+    fvc::makeRelative(phiHbyA, U);
+    adjustPhi(phiHbyA, U, p);
+    fvc::makeAbsolute(phiHbyA, U);
+}
+
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, U, phiHbyA, rAU, MRF);
+
+// Non-orthogonal pressure corrector loop
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rAU, p)
+     ==
+        fvc::div(phiHbyA)
+    );
+
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA - pEqn.flux();
+    }
+}
+
+#include "continuityErrs.H"
+
+p.relax();
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+
+// Correct rhoUf if the mesh is moving
+fvc::correctUf(Uf, U, phi);
+
+// Make the fluxes relative to the mesh motion
+fvc::makeRelative(phi, U);
+
--- a/applications/solvers/lagrangian/parcelFoam/parcelFoam.C
+++ b/applications/solvers/lagrangian/parcelFoam/parcelFoam.C
@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    parcelFoam
+
+Group
+    grpLagrangianSolvers
+
+Description
+    Transient solver for incompressible, turbulent flow with kinematic,
+    particle cloud, and surface film modelling.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "surfaceFilmModel.H"
+#include "parcelCloudModelList.H"
+#include "fvOptions.H"
+#include "pimpleControl.H"
+#include "CorrectPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Transient solver for incompressible, turbulent flow"
+        " with kinematic particle clouds"
+        " and surface film modelling."
+    );
+
+    #define CREATE_MESH createMeshesPostProcess.H
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createUfIfPresent.H"
+
+    turbulence->validate();
+
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readDyMControls.H"
+        #include "CourantNo.H"
+        #include "setMultiRegionDeltaT.H"
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Store the particle positions
+        parcels.storeGlobalPositions();
+
+        // Do any mesh changes
+        mesh.update();
+
+        if (pimple.solveFlow() && mesh.changing())
+        {
+            MRF.update();
+
+            if (correctPhi)
+            {
+                // Calculate absolute flux
+                // from the mapped surface velocity
+                phi = mesh.Sf() & Uf();
+
+                #include "../../incompressible/pimpleFoam/correctPhi.H"
+
+                // Make the fluxes relative to the mesh-motion
+                fvc::makeRelative(phi, U);
+            }
+
+            if (checkMeshCourantNo)
+            {
+                #include "meshCourantNo.H"
+            }
+        }
+
+        parcels.evolve();
+        surfaceFilm.evolve();
+
+        if (pimple.solveFlow())
+        {
+            // --- PIMPLE loop
+            while (pimple.loop())
+            {
+                #include "UEqn.H"
+
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    laminarTransport.correct();
+                    turbulence->correct();
+                }
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingTypeCloud parcels
+reactingTypeCloud parcels
 (
    "reactingCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/reactingHeterogenousParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/reactingHeterogenousParcelFoam.C
@ -37,8 +37,8 @@ Description

 \*---------------------------------------------------------------------------*/

-#define CLOUD_BASE_TYPE HeterogeneousReacting
-#define CLOUD_BASE_TYPE_NAME "HeterogeneousReacting"
+#define CLOUD_BASE_TYPE heterogeneousReacting
+#define CLOUD_BASE_TYPE_NAME "heterogeneousReacting"

 #include "reactingParcelFoam.C"

--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@ -53,12 +53,12 @@ Description
 #include "cloudMacros.H"

 #ifndef CLOUD_BASE_TYPE
-    #define CLOUD_BASE_TYPE ReactingMultiphase
+    #define CLOUD_BASE_TYPE reactingMultiphase
    #define CLOUD_BASE_TYPE_NAME "reacting"
 #endif

 #include CLOUD_INCLUDE_FILE(CLOUD_BASE_TYPE)
-#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
+#define reactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -81,7 +81,6 @@ int main(int argc, char *argv[])
    #include "createDyMControls.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
-    #include "createRegionControls.H"
    #include "initContinuityErrs.H"
    #include "createRhoUfIfPresent.H"

@ -105,7 +104,7 @@ int main(int argc, char *argv[])
        // so that it can be mapped and used in correctPhi
        // to ensure the corrected phi has the same divergence
        autoPtr<volScalarField> divrhoU;
-        if (solvePrimaryRegion && correctPhi)
+        if (pimple.solveFlow() && correctPhi)
        {
            divrhoU.reset
            (
@ -133,7 +132,7 @@ int main(int argc, char *argv[])

        // Store momentum to set rhoUf for introduced faces.
        autoPtr<volVectorField> rhoU;
-        if (solvePrimaryRegion && rhoUf.valid())
+        if (pimple.solveFlow() && rhoUf.valid())
        {
            rhoU.reset(new volVectorField("rhoU", rho*U));
        }
@ -144,7 +143,7 @@ int main(int argc, char *argv[])
        // Do any mesh changes
        mesh.update();

-        if (solvePrimaryRegion && mesh.changing())
+        if (pimple.solveFlow() && mesh.changing())
        {
            gh = (g & mesh.C()) - ghRef;
            ghf = (g & mesh.Cf()) - ghRef;
@ -172,7 +171,7 @@ int main(int argc, char *argv[])
        parcels.evolve();
        surfaceFilm.evolve();

-        if (solvePrimaryRegion)
+        if (pimple.solveFlow())
        {
            if (pimple.nCorrPIMPLE() <= 1)
            {
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
@ -7,7 +7,6 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing " << CLOUD_BASE_TYPE_NAME << " cloud" << endl;
-basicReactingTypeCloud parcels
+reactingTypeCloud parcels
 (
    word(CLOUD_BASE_TYPE_NAME) & "Cloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
@ -45,13 +45,13 @@ Description
 #include "cloudMacros.H"

 #ifndef CLOUD_BASE_TYPE
-    #define CLOUD_BASE_TYPE ReactingMultiphase
+    #define CLOUD_BASE_TYPE reactingMultiphase
    //#define CLOUD_BASE_TYPE_NAME "reactingMultiphase" Backwards compat
    #define CLOUD_BASE_TYPE_NAME "reacting"
 #endif

 #include CLOUD_INCLUDE_FILE(CLOUD_BASE_TYPE)
-#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
+#define reactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/lagrangian/rhoParcelFoam/EEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/EEqn.H
@ -0,0 +1,41 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        rho*(U&g)
+      + parcels.Sh(he)
+      + surfaceFilm.Sh()
+      + radiation->Sh(thermo, he)
+      + Qdot
+      + fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
--- a/applications/solvers/lagrangian/rhoParcelFoam/Make/files
+++ b/applications/solvers/lagrangian/rhoParcelFoam/Make/files
@ -0,0 +1,3 @@
+rhoParcelFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/rhoParcelFoam
--- a/applications/solvers/lagrangian/rhoParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/rhoParcelFoam/Make/options
@ -1,52 +1,58 @@
 EXE_INC = \
+    -I../reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-    -I$(LIB_SRC)/lagrangian/turbulence/lnInclude \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/faOptions/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam

-LIB_LIBS = \
+EXE_LIBS = \
    -lfiniteVolume \
-    -lsurfMesh \
+    -lfvOptions \
    -lmeshTools \
-    -llagrangian \
-    -llagrangianIntermediate \
-    -llagrangianTurbulence \
-    -ldistributionModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -lsampling \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lthermoTools \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
-    -lthermophysicalProperties \
    -lreactionThermophysicalModels \
+    -lthermophysicalProperties \
    -lSLGThermo \
-    -lradiationModels \
-    -lturbulenceModels \
-    -lincompressibleTurbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lincompressibleTransportModels \
+    -lchemistryModel \
    -lregionModels \
+    -lradiationModels \
    -lsurfaceFilmModels \
-    -ldynamicMesh \
-    -ldynamicFvMesh \
+    -lsurfaceFilmDerivedFvPatchFields \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lODE \
+    -lcombustionModels \
    -lregionFaModels \
-    -lfiniteArea
+    -lfiniteArea \
+    -lfaOptions
--- a/applications/solvers/lagrangian/rhoParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/UEqn.H
@ -0,0 +1,34 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        parcels.SU(U)
+      + fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+          ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
+    }
--- a/applications/solvers/lagrangian/rhoParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/YEqn.H
@ -0,0 +1,51 @@
+tmp<fv::convectionScheme<scalar>> mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+
+{
+    combustion->correct();
+    Qdot = combustion->Qdot();
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (i != inertIndex && composition.active(i))
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YEqn
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+              ==
+                parcels.SYi(i, Yi)
+              + fvOptions(rho, Yi)
+              + combustion->R(Yi)
+              + surfaceFilm.Srho(i)
+            );
+
+            YEqn.relax();
+
+            fvOptions.constrain(YEqn);
+
+            YEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}
--- a/applications/solvers/lagrangian/rhoParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createClouds.H
@ -0,0 +1,3 @@
+Info<< "\nConstructing cloud" << endl;
+
+auto parcels = parcelCloudModelList(g, rho, U, slgThermo);
--- a/applications/solvers/lagrangian/rhoParcelFoam/createFieldRefs.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createFieldRefs.H
@ -0,0 +1,5 @@
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
+const label inertIndex(composition.species().find(inertSpecie));
+
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
--- a/applications/solvers/lagrangian/rhoParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createFields.H
@ -0,0 +1,132 @@
+#include "createRDeltaT.H"
+
+Info<< "Reading thermophysical properties\n" << endl;
+autoPtr<rhoReactionThermo> pThermo(rhoReactionThermo::New(mesh));
+rhoReactionThermo& thermo = pThermo();
+thermo.validate(args.executable(), "h", "e");
+
+SLGThermo slgThermo(mesh, thermo);
+
+basicSpecieMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+const word inertSpecie(thermo.get<word>("inertSpecie"));
+if
+(
+    !composition.species().found(inertSpecie)
+  && composition.species().size() > 0
+)
+{
+    FatalIOErrorIn(args.executable().c_str(), thermo)
+        << "Inert specie " << inertSpecie << " not found in available species "
+        << composition.species()
+        << exit(FatalIOError);
+}
+
+Info<< "Creating field rho\n" << endl;
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+volScalarField& p = thermo.p();
+
+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 combustion model\n" << endl;
+autoPtr<CombustionModel<rhoReactionThermo>> combustion
+(
+    CombustionModel<rhoReactionThermo>::New(thermo, turbulence())
+);
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Force p_rgh to be consistent with p
+p_rgh = p - rho*gh;
+
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
+mesh.setFluxRequired(p_rgh.name());
+
+Info<< "Creating multi-variate interpolation scheme\n" << endl;
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField Qdot
+(
+    IOobject
+    (
+        "Qdot",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimEnergy/dimVolume/dimTime, Zero)
+);
+
+#include "createDpdt.H"
+
+#include "createK.H"
+
+#include "createMRF.H"
+#include "createRadiationModel.H"
+#include "createClouds.H"
+#include "createSurfaceFilmModel.H"
+#include "createFvOptions.H"
--- a/applications/solvers/lagrangian/rhoParcelFoam/createMeshesPostProcess.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createMeshesPostProcess.H
@ -0,0 +1,35 @@
+#include "createMesh.H"
+
+dictionary filmDict;
+
+IOobject io
+(
+    "surfaceFilmProperties",
+    mesh.time().constant(),
+    mesh,
+    IOobject::READ_IF_PRESENT,
+    IOobject::NO_WRITE,
+    false
+);
+
+if (io.typeHeaderOk<IOdictionary>())
+{
+    IOdictionary propDict(io);
+
+    filmDict = std::move(propDict);
+
+    const word filmRegionName = filmDict.get<word>("region");
+
+    fvMesh filmMesh
+    (
+        IOobject
+        (
+            filmRegionName,
+            runTime.timeName(),
+            runTime,
+            IOobject::MUST_READ
+        )
+    );
+}
+
+
--- a/applications/solvers/lagrangian/rhoParcelFoam/createSurfaceFilmModel.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/createSurfaceFilmModel.H
@ -0,0 +1,6 @@
+Info<< "\nConstructing surface film model" << endl;
+
+autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
+(
+    regionModels::surfaceFilmModel::New(mesh, g)
+);
--- a/applications/solvers/lagrangian/rhoParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/pEqn.H
@ -0,0 +1,96 @@
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
+    )
+  + phig
+);
+
+fvc::makeRelative(phiHbyA, rho, U);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+  + fvc::div(phiHbyA)
+ ==
+    parcels.Srho()
+  + surfaceFilm.Srho()
+  + fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix p_rghEqn
+    (
+        p_rghDDtEqn
+      - fvm::laplacian(rhorAUf, p_rgh)
+    );
+
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA + p_rghEqn.flux();
+
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
+
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
+    }
+}
+
+p = p_rgh + rho*gh;
+
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+    p_rgh = p - rho*gh;
+}
+else if (pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Correct rhoUf if the mesh is moving
+fvc::correctRhoUf(rhoUf, rho, U, phi);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+
+    if (mesh.moving())
+    {
+        dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
+    }
+}
--- a/applications/solvers/lagrangian/rhoParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/rhoEqn.H
@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        parcels.Srho(rho)
+      + surfaceFilm.Srho()
+      + fvOptions(rho)
+    );
+
+    rhoEqn.solve();
+
+    fvOptions.correct(rho);
+}
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/rhoParcelFoam/rhoParcelFoam.C
+++ b/applications/solvers/lagrangian/rhoParcelFoam/rhoParcelFoam.C
@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2020 OpenFOAM Foundation
+    Copyright (C) 2018-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rhoParcelFoam
+
+Group
+    grpLagrangianSolvers
+
+Description
+    Transient solver for compressible, turbulent flow with a reacting,
+    multiphase particle cloud, and surface film modelling.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "turbulentFluidThermoModel.H"
+#include "surfaceFilmModel.H"
+#include "rhoReactionThermo.H"
+#include "CombustionModel.H"
+#include "radiationModel.H"
+#include "SLGThermo.H"
+#include "fvOptions.H"
+#include "pimpleControl.H"
+#include "pressureControl.H"
+#include "CorrectPhi.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+#include "parcelCloudModelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Transient solver for compressible, turbulent flow"
+        " with lagrangian parcel clouds and surface film modelling."
+    );
+
+    #define CREATE_MESH createMeshesPostProcess.H
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createDyMControls.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "initContinuityErrs.H"
+    #include "createRhoUfIfPresent.H"
+
+    turbulence->validate();
+
+    if (!LTS)
+    {
+        #include "compressibleCourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readDyMControls.H"
+
+        // Store divrhoU from the previous mesh
+        // so that it can be mapped and used in correctPhi
+        // to ensure the corrected phi has the same divergence
+        autoPtr<volScalarField> divrhoU;
+        if (pimple.solveFlow() && correctPhi)
+        {
+            divrhoU.reset
+            (
+                new volScalarField
+                (
+                    "divrhoU",
+                    fvc::div(fvc::absolute(phi, rho, U))
+                )
+            );
+        }
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "compressibleCourantNo.H"
+            #include "setMultiRegionDeltaT.H"
+        }
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Store momentum to set rhoUf for introduced faces.
+        autoPtr<volVectorField> rhoU;
+        if (pimple.solveFlow() && rhoUf.valid())
+        {
+            rhoU.reset(new volVectorField("rhoU", rho*U));
+        }
+
+        // Store the particle positions
+        parcels.storeGlobalPositions();
+
+        // Do any mesh changes
+        mesh.update();
+
+        if (pimple.solveFlow() && mesh.changing())
+        {
+            gh = (g & mesh.C()) - ghRef;
+            ghf = (g & mesh.Cf()) - ghRef;
+
+            MRF.update();
+
+            if (correctPhi)
+            {
+                // Calculate absolute flux
+                // from the mapped surface velocity
+                phi = mesh.Sf() & rhoUf();
+
+                #include "../../compressible/rhoPimpleFoam/correctPhi.H"
+
+                // Make the fluxes relative to the mesh-motion
+                fvc::makeRelative(phi, rho, U);
+            }
+
+            if (checkMeshCourantNo)
+            {
+                #include "meshCourantNo.H"
+            }
+        }
+
+        parcels.evolve();
+        surfaceFilm.evolve();
+
+        if (pimple.solveFlow())
+        {
+            if (pimple.nCorrPIMPLE() <= 1)
+            {
+                #include "rhoEqn.H"
+            }
+
+            // --- PIMPLE loop
+            while (pimple.loop())
+            {
+                #include "UEqn.H"
+                #include "YEqn.H"
+                #include "EEqn.H"
+
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
+
+            rho = thermo.rho();
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/rhoParcelFoam/setMultiRegionDeltaT.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/setMultiRegionDeltaT.H
@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    setMultiRegionDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum Courant numbers.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    const scalar maxDeltaTFact =
+        min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
+    const scalar deltaTFact =
+        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            deltaTFact*runTime.deltaTValue(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/rhoParcelFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/setRDeltaT.H
@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    volScalarField& rDeltaT = trDeltaT.ref();
+
+    const dictionary& pimpleDict = pimple.dict();
+
+    // Maximum flow Courant number
+    scalar maxCo(pimpleDict.get<scalar>("maxCo"));
+
+    // Maximum time scale
+    scalar maxDeltaT(pimpleDict.getOrDefault<scalar>("maxDeltaT", GREAT));
+
+    // Smoothing parameter (0-1) when smoothing iterations > 0
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.getOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    // Damping coefficient (1-0)
+    scalar rDeltaTDampingCoeff
+    (
+        pimpleDict.getOrDefault<scalar>("rDeltaTDampingCoeff", 0.2)
+    );
+
+    // Maximum change in cell temperature per iteration
+    // (relative to previous value)
+    scalar alphaTemp(pimpleDict.getOrDefault("alphaTemp", 0.05));
+
+
+    Info<< "Time scales min/max:" << endl;
+
+    // Cache old reciprocal time scale field
+    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Flow time scale
+    {
+        rDeltaT.ref() =
+        (
+            fvc::surfaceSum(mag(phi))()()
+           /((2*maxCo)*mesh.V()*rho())
+        );
+
+        // Limit the largest time scale
+        rDeltaT.max(1/maxDeltaT);
+
+        Info<< "    Flow        = "
+            << gMin(1/rDeltaT.primitiveField()) << ", "
+            << gMax(1/rDeltaT.primitiveField()) << endl;
+    }
+
+    // Reaction source time scale
+    {
+        volScalarField::Internal rDeltaTT
+        (
+            mag
+            (
+                parcels.hsTrans()/(mesh.V()*runTime.deltaT())
+              + Qdot()
+            )
+           /(
+               alphaTemp
+              *rho()
+              *thermo.Cp()()()
+              *T()
+           )
+        );
+
+        Info<< "    Temperature = "
+            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
+            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+
+        rDeltaT.ref() = max
+        (
+            rDeltaT(),
+            rDeltaTT
+        );
+    }
+
+    // Update the boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    // Spatially smooth the time scale field
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+    }
+
+    Info<< "    Overall     = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
@ -7,7 +7,6 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@ -47,7 +46,6 @@ EXE_LIBS = \
    -lsurfaceFilmModels \
    -lcombustionModels \
    -lsampling \
-    -lcoalCombustion \
    -lregionFaModels \
    -lfiniteArea \
    -lfaOptions
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing coal cloud" << endl;
-coalCloud parcels
+reactingMultiphaseCloud parcels
 (
    "reactingCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/simpleCoalParcelFoam.C
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/simpleCoalParcelFoam.C
@ -37,7 +37,7 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "coalCloud.H"
+#include "reactingMultiphaseCloud.H"
 #include "rhoReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/sprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/Make/options
@ -35,7 +35,6 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -llagrangianSpray \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/lagrangian/sprayFoam/createClouds.H
+++ b/applications/solvers/lagrangian/sprayFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicSprayCloud parcels
+sprayCloud parcels
 (
    "sprayCloud",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
@ -39,7 +39,6 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -llagrangianSpray \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/engineFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/engineFoam.C
@ -37,7 +37,7 @@ Description
 #include "engineTime.H"
 #include "engineMesh.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicSprayCloud.H"
+#include "sprayCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
@ -34,7 +34,6 @@ EXE_LIBS = \
    -lthermoTools \
    -llagrangian \
    -llagrangianIntermediate \
-    -llagrangianSpray \
    -llagrangianTurbulence \
    -lspecie \
    -lcompressibleTransportModels \
--- a/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/simpleSprayFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/simpleSprayFoam.C
@ -35,7 +35,7 @@ Description

 \*---------------------------------------------------------------------------*/

-#define CLOUD_BASE_TYPE Spray
+#define CLOUD_BASE_TYPE spray
 #define CLOUD_BASE_TYPE_NAME "spray"

 #include "simpleReactingParcelFoam.C"
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
@ -40,7 +40,6 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -llagrangianSpray \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C
@ -38,7 +38,7 @@ Description
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
 #include "turbulenceModel.H"
-#include "basicSprayCloud.H"
+#include "sprayCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayFoam.C
@ -37,7 +37,7 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicSprayCloud.H"
+#include "sprayCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H
@ -55,11 +55,11 @@ const word kinematicCloudName
 );

 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-basicKinematicCloud kinematicCloud
+kinematicCloud kCloud
 (
    kinematicCloudName,
+    g,
    rho,
    U,
-    thermo.mu(),
-    g
+    thermo.mu()
 );
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C
@ -37,7 +37,7 @@ Description
 #include "dynamicFvMesh.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicKinematicCloud.H"
+#include "kinematicCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -73,14 +73,14 @@ int main(int argc, char *argv[])
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

-        kinematicCloud.storeGlobalPositions();
+        kCloud.storeGlobalPositions();

        mesh.update();

        U.correctBoundaryConditions();

-        Info<< "Evolving " << kinematicCloud.name() << endl;
-        kinematicCloud.evolve();
+        Info<< "Evolving " << kCloud.name() << endl;
+        kCloud.evolve();

        runTime.write();

--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
@ -39,7 +39,7 @@ Description
 #include "fvCFD.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicKinematicCloud.H"
+#include "kinematicCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -75,8 +75,8 @@ int main(int argc, char *argv[])
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

-        Info<< "Evolving " << kinematicCloud.name() << endl;
-        kinematicCloud.evolve();
+        Info<< "Evolving " << kCloud.name() << endl;
+        kCloud.evolve();

        runTime.write();

--- a/applications/solvers/multiphase/MPPICInterFoam/MPPICInterFoam.C
+++ b/applications/solvers/multiphase/MPPICInterFoam/MPPICInterFoam.C
@ -52,7 +52,7 @@ Description
 #include "CorrectPhi.H"
 #include "fvcSmooth.H"

-#include "basicKinematicCloud.H"
+#include "kinematicCloud.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -111,12 +111,12 @@ int main(int argc, char *argv[])

        Info<< "Time = " << runTime.timeName() << nl << endl;

-        Info<< "Evolving " << kinematicCloud.name() << endl;
+        Info<< "Evolving " << kCloud.name() << endl;

-        kinematicCloud.evolve();
+        kCloud.evolve();

        // Update continuous phase volume fraction field
-        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+        alphac = max(1.0 - kCloud.theta(), alphacMin);
        alphac.correctBoundaryConditions();

         Info<< "Continuous phase-1 volume fraction = "
@ -130,7 +130,7 @@ int main(int argc, char *argv[])
        alphaPhic = alphacf*phi;
        alphacRho = alphac*rho;

-        fvVectorMatrix cloudSU(kinematicCloud.SU(U));
+        fvVectorMatrix cloudSU(kCloud.SU(U));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/multiphase/MPPICInterFoam/Make/options
+++ b/applications/solvers/multiphase/MPPICInterFoam/Make/options
@ -12,11 +12,14 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
--- a/applications/solvers/multiphase/MPPICInterFoam/createFields.H
+++ b/applications/solvers/multiphase/MPPICInterFoam/createFields.H
@ -147,13 +147,13 @@ volScalarField alphacRho(alphac*rho);
 alphacRho.oldTime();

 Info<< "Constructing kinematicCloud " << endl;
-basicKinematicCloud kinematicCloud
+kinematicCloud kCloud
 (
    "kinematicCloud",
+    g,
    rho,
    U,
-    mu,
-    g
+    mu
 );

 // Particle fraction upper limit
@ -161,13 +161,13 @@ scalar alphacMin
 (
    1.0
  - (
-        kinematicCloud.particleProperties().subDict("constantProperties")
+        kCloud.particleProperties().subDict("constantProperties")
       .get<scalar>("alphaMax")
    )
 );

 // Update alphac from the particle locations
-alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+alphac = max(1.0 - kCloud.theta(), alphacMin);
 alphac.correctBoundaryConditions();

 surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3DKernel.H
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3DKernel.H
@ -30,8 +30,8 @@ Description

 \*---------------------------------------------------------------------------*/

-#ifndef Foam_CGALTriangulation3DKernel_H
-#define Foam_CGALTriangulation3DKernel_H
+#ifndef CGALTriangulation3DKernel_H
+#define CGALTriangulation3DKernel_H

 // Silence boost bind deprecation warnings (before CGAL-5.2.1)
 #include "CGAL/version.h"
@ -54,19 +54,9 @@ Description
 //    #include "CGAL/Robust_circumcenter_traits_3.h"
 //    typedef CGAL::Robust_circumcenter_traits_3<baseK>           K;

-#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050500000)
-    // Prior to CGAL-5.5
-
    #include "CGAL/Robust_circumcenter_filtered_traits_3.h"
    typedef CGAL::Robust_circumcenter_filtered_traits_3<baseK>  K;

-#else
-
-    #include "CGAL/Robust_weighted_circumcenter_filtered_traits_3.h"
-    typedef CGAL::Robust_weighted_circumcenter_filtered_traits_3<baseK> K;
-
-#endif
-
 #else

    // Very robust but expensive kernel
--- a/applications/utilities/preProcessing/viewFactorsGen/Make/options
+++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options
@ -4,6 +4,7 @@ EXE_INC = \
    -Wno-old-style-cast \
    $(COMP_FLAGS) \
    ${CGAL_INC} \
+    -DCGAL_HEADER_ONLY \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/surfMesh/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
@ -11,6 +12,7 @@ EXE_INC = \


 EXE_LIBS = \
+    /* ${CGAL_LIBS} */ \
    -lfiniteVolume \
    -lsurfMesh \
    -lmeshTools \
--- a/etc/bashrc
+++ b/etc/bashrc
@ -55,7 +55,7 @@

 # [WM_PROJECT_VERSION] - A human-readable version name
 # A development version is often named 'com' - as in www.openfoam.com
-export WM_PROJECT_VERSION=v2212
+export WM_PROJECT_VERSION=com

 #------------------------------------------------------------------------------
 # Configuration environment variables.
--- a/etc/caseDicts/postProcessing/solvers/scalarTransport/scalarTransport
+++ b/etc/caseDicts/postProcessing/solvers/scalarTransport/scalarTransport
@ -18,8 +18,6 @@ Description
 type            scalarTransport;
 libs            ("libsolverFunctionObjects.so");

-writeControl    writeTime;
-
 field           s;
 schemesField    s;
 D               1e-09;
--- a/etc/cshrc
+++ b/etc/cshrc
@ -55,7 +55,7 @@

 # [WM_PROJECT_VERSION] - A human-readable version name
 # A development version is often named 'com' - as in www.openfoam.com
-setenv WM_PROJECT_VERSION v2212
+setenv WM_PROJECT_VERSION com

 #------------------------------------------------------------------------------
 # Configuration environment variables.
--- a/src/OpenFOAM/db/dictionary/functionEntries/evalEntry/evalEntry.C
+++ b/src/OpenFOAM/db/dictionary/functionEntries/evalEntry/evalEntry.C
@ -174,7 +174,7 @@ Foam::tokenList Foam::functionEntries::evalEntry::evaluate
        result.writeField(toks);
    }

-    return tokenList(std::move(toks.tokens()));
+    return std::move(toks);
 }


--- a/src/OpenFOAM/expressions/exprDriver/exprDriverFields.C
+++ b/src/OpenFOAM/expressions/exprDriver/exprDriverFields.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -41,11 +41,9 @@ void Foam::expressions::exprDriver::fill_random
 {
    if (seed <= 0)
    {
-        const TimeState* ts = this->timeState();
-
-        if (ts)
+        if (timeStatePtr_)
        {
-            seed = (ts->timeIndex() - seed);
+            seed = (timeStatePtr_->timeIndex() - seed);
        }
        else
        {
--- a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
@ -152,11 +152,10 @@ slicedBoundaryField
                new SlicedPatchField<Type>
                (
                    mesh.boundary()[patchi],
-                    DimensionedField<Type, GeoMesh>::null()
+                    DimensionedField<Type, GeoMesh>::null(),
+                    bField[patchi]
                )
            );
-
-            bf[patchi].UList<Type>::shallowCopy(bField[patchi]);
        }
    }

--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@ -201,8 +201,7 @@ inline Foam::Matrix<Form, Type>::Matrix
 )
 :
    mRows_(Mb.m()),
-    nCols_(Mb.n()),
-    v_(nullptr)
+    nCols_(Mb.n())
 {
    doAlloc();

@ -224,8 +223,7 @@ inline Foam::Matrix<Form, Type>::Matrix
 )
 :
    mRows_(Mb.m()),
-    nCols_(Mb.n()),
-    v_(nullptr)
+    nCols_(Mb.n())
 {
    doAlloc();

--- a/src/OpenFOAM/parallel/globalIndex/globalIndex.H
+++ b/src/OpenFOAM/parallel/globalIndex/globalIndex.H
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -489,12 +489,12 @@ public:

        //- Collect indirect data in processor order on master
        //  Handles contiguous/non-contiguous data, skips empty fields.
-        template<class ProcIDsContainer, class Type, class Addr>
+        template<class Type, class Addr>
        static void gather
        (
            const labelUList& offsets,  //!< offsets (master only)
            const label comm,           //!< communicator
-            const ProcIDsContainer& procIDs,
+            const UList<int>& procIDs,
            const IndirectListBase<Type, Addr>& fld,
            List<Type>& allFld,         //! output field (master only)
            const int tag = UPstream::msgType(),
--- a/src/OpenFOAM/parallel/globalIndex/globalIndexTemplates.C
+++ b/src/OpenFOAM/parallel/globalIndex/globalIndexTemplates.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2013-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -248,12 +248,12 @@ void Foam::globalIndex::gather
 }


-template<class ProcIDsContainer, class Type, class Addr>
+template<class Type, class Addr>
 void Foam::globalIndex::gather
 (
    const labelUList& off,  // needed on master only
    const label comm,
-    const ProcIDsContainer& procIDs,
+    const UList<int>& procIDs,
    const IndirectListBase<Type, Addr>& fld,
    List<Type>& allFld,
    const int tag,
@ -368,7 +368,7 @@ void Foam::globalIndex::gather
        (
            offsets_,  // needed on master only
            comm,
-            UPstream::allProcs(comm),   // All communicator ranks
+            UPstream::procID(comm),
            sendData,
            allData,
            tag,
@ -404,7 +404,7 @@ void Foam::globalIndex::gather
        (
            offsets_,  // needed on master only
            comm,
-            UPstream::allProcs(comm),   // All communicator ranks
+            UPstream::procID(comm),
            sendData,
            allData,
            tag,
@ -622,7 +622,7 @@ void Foam::globalIndex::mpiGather
            (
                offsets_,  // needed on master only
                comm,
-                UPstream::allProcs(comm),   // All communicator ranks
+                UPstream::procID(comm),
                sendData,
                allData,
                tag,
@ -967,7 +967,7 @@ void Foam::globalIndex::scatter
        (
            offsets_,  // needed on master only
            comm,
-            UPstream::allProcs(comm),  // All communicator ranks
+            UPstream::procID(comm),
            allData,
            localData,
            tag,
--- a/src/OpenFOAM/primitives/Tensor/TensorI.H
+++ b/src/OpenFOAM/primitives/Tensor/TensorI.H
@ -409,10 +409,6 @@ inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::T() const


 template<class Cmpt>
-#if defined(__GNUC__) && !defined(__clang__)
-// Workaround for gcc (11+) that fails to handle tensor dot vector
-__attribute__((optimize("no-tree-vectorize")))
-#endif
 inline Foam::Tensor<Cmpt>
 Foam::Tensor<Cmpt>::inner(const Tensor<Cmpt>& t2) const
 {
@ -436,10 +432,6 @@ Foam::Tensor<Cmpt>::inner(const Tensor<Cmpt>& t2) const


 template<class Cmpt>
-#if defined(__GNUC__) && !defined(__clang__)
-// Workaround for gcc (11+) that fails to handle tensor dot vector
-__attribute__((optimize("no-tree-vectorize")))
-#endif
 inline Foam::Tensor<Cmpt>
 Foam::Tensor<Cmpt>::schur(const Tensor<Cmpt>& t2) const
 {
@ -875,10 +867,6 @@ operator&(const Tensor<Cmpt>& t1, const Tensor<Cmpt>& t2)

 //- Inner-product of a SphericalTensor and a Tensor
 template<class Cmpt>
-#if defined(__GNUC__) && !defined(__clang__)
-// Workaround for gcc (11+) that fails to handle tensor dot vector
-__attribute__((optimize("no-tree-vectorize")))
-#endif
 inline Tensor<Cmpt>
 operator&(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 {
@ -893,10 +881,6 @@ operator&(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)

 //- Inner-product of a Tensor and a SphericalTensor
 template<class Cmpt>
-#if defined(__GNUC__) && !defined(__clang__)
-// Workaround for gcc (11+) that fails to handle tensor dot vector
-__attribute__((optimize("no-tree-vectorize")))
-#endif
 inline Tensor<Cmpt>
 operator&(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
 {
@ -911,10 +895,6 @@ operator&(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)

 //- Inner-product of a SymmTensor and a Tensor
 template<class Cmpt>
-#if defined(__GNUC__) && !defined(__clang__)
-// Workaround for gcc (11+) that fails to handle tensor dot vector
-__attribute__((optimize("no-tree-vectorize")))
-#endif
 inline Tensor<Cmpt>
 operator&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 {
@ -937,10 +917,6 @@ operator&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)

 //- Inner-product of a Tensor and a SymmTensor
 template<class Cmpt>
-#if defined(__GNUC__) && !defined(__clang__)
-// Workaround for gcc (11+) that fails to handle tensor dot vector
-__attribute__((optimize("no-tree-vectorize")))
-#endif
 inline Tensor<Cmpt>
 operator&(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
 {
@ -981,10 +957,6 @@ operator&(const Tensor<Cmpt>& t, const Vector<Cmpt>& v)

 //- Inner-product of a Vector and a Tensor
 template<class Cmpt>
-#if defined(__GNUC__) && !defined(__clang__)
-// Workaround for gcc (11+) that fails to handle tensor dot vector
-__attribute__((optimize("no-tree-vectorize")))
-#endif
 inline Vector<Cmpt>
 operator&(const Vector<Cmpt>& v, const Tensor<Cmpt>& t)
 {
--- a/src/OpenFOAM/primitives/strings/word/word.C
+++ b/src/OpenFOAM/primitives/strings/word/word.C
@ -29,7 +29,6 @@ License
 #include "word.H"
 #include "debug.H"
 #include <cctype>
-#include <cstdint>
 #include <sstream>

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
--- a/src/Pstream/mpi/UPstream.C
+++ b/src/Pstream/mpi/UPstream.C
@ -641,16 +641,7 @@ void Foam::UPstream::freePstreamCommunicator(const label communicator)
    }

    // Not touching the first two communicators (SELF, WORLD)
-    // or anything out-of bounds.
-    //
-    // No UPstream communicator indices when MPI is initialized outside
-    // of OpenFOAM - thus needs a bounds check too!
-
-    if
-    (
-        communicator > 1
-     && (communicator < PstreamGlobals::MPICommunicators_.size())
-    )
+    if (communicator > 1)
    {
        if (MPI_COMM_NULL != PstreamGlobals::MPICommunicators_[communicator])
        {
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2022-2023 OpenCFD Ltd.
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -131,9 +131,6 @@ makeLESModel(dynamicKEqn);
 #include "dynamicLagrangian.H"
 makeLESModel(dynamicLagrangian);

-#include "sigma.H"
-makeLESModel(sigma);
-
 #include "SpalartAllmarasDES.H"
 makeLESModel(SpalartAllmarasDES);

--- a/src/TurbulenceModels/turbulenceModels/laminar/generalizedNewtonian/generalizedNewtonianViscosityModels/Casson/Casson.C
+++ b/src/TurbulenceModels/turbulenceModels/laminar/generalizedNewtonian/generalizedNewtonianViscosityModels/Casson/Casson.C
@ -79,8 +79,8 @@ bool Foam::laminarModels::generalizedNewtonianViscosityModels::Casson::read

    coeffs.readEntry("m", m_);
    coeffs.readEntry("tau0", tau0_);
-    coeffs.readEntry("nuMin", nuMin_);
-    coeffs.readEntry("nuMax", nuMax_);
+    coeffs.readEntry("nuMin_", nuMin_);
+    coeffs.readEntry("nuMax_", nuMax_);

    return true;
 }
--- a/src/fileFormats/nastran/NASCore.C
+++ b/src/fileFormats/nastran/NASCore.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2024 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -122,8 +122,7 @@ std::string Foam::fileFormats::NASCore::nextNasField
 (
    const std::string& str,
    std::string::size_type& pos,
-    const std::string::size_type width,
-    const bool free_format
+    std::string::size_type len
 )
 {
    const auto beg = pos;
@ -131,23 +130,15 @@ std::string Foam::fileFormats::NASCore::nextNasField

    if (end == std::string::npos)
    {
-        if (free_format)
-        {
-            // Nothing left
-            pos = str.size();
-            return str.substr(beg);
-        }
-
-        // Fixed format - continue after field width
-        pos = beg + width;
-        return str.substr(beg, width);
+        pos = beg + len;    // Continue after field width
    }
    else
    {
-        // Free format - continue after comma
-        pos = end + 1;
-        return str.substr(beg, (end - beg));
+        len = (end - beg);  // Efffective width
+        pos = end + 1;      // Continue after comma
    }
+
+    return str.substr(beg, len);
 }


@ -244,8 +235,8 @@ void Foam::fileFormats::NASCore::writeCoord
    // 2 ID   : point ID - requires starting index of 1
    // 3 CP   : coordinate system ID                (blank)
    // 4 X1   : point x coordinate
-    // 5 X2   : point y coordinate
-    // 6 X3   : point z coordinate
+    // 5 X2   : point x coordinate
+    // 6 X3   : point x coordinate
    // 7 CD   : coordinate system for displacements (blank)
    // 8 PS   : single point constraints            (blank)
    // 9 SEID : super-element ID
--- a/src/fileFormats/nastran/NASCore.H
+++ b/src/fileFormats/nastran/NASCore.H
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2024 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -48,6 +48,7 @@ SourceFiles

 namespace Foam
 {
+
 namespace fileFormats
 {

@ -73,18 +74,18 @@ public:
        //- Output load format
        enum loadFormat
        {
-            PLOAD2,     //!< Face load (eg, pressure)
-            PLOAD4      //!< Vertex load
+            PLOAD2,
+            PLOAD4
        };

-        //- Selection names for the NASTRAN load formats
+        //- Selection names for the NASTRAN file field formats
        static const Enum<loadFormat> loadFormatNames;


    // Constructors

        //- Default construct
-        NASCore() noexcept = default;
+        NASCore() = default;


    // Public Static Member Functions
@ -92,20 +93,18 @@ public:
        //- Extract numbers from things like "-2.358-8" (same as "-2.358e-8")
        static scalar readNasScalar(const std::string& str);

-        //- A std::string::substr() variant to handle fixed-format and
-        //- free-format NASTRAN.
-        //  Returns the substr until the next comma (if found)
-        //  or the given fixed width
+        //- A string::substr() to handle fixed-format and free-format NASTRAN.
+        //  Returns the substr to the next comma (if found) or the given length
+        //
+        //  \param str The string to extract from
+        //  \param pos On input, the position of the first character of the
+        //      substring. On output, advances to the next position to use.
+        //  \param len The fixed-format length to use if a comma is not found.
        static std::string nextNasField
        (
-            //! The string to extract from
            const std::string& str,
-            //! [in,out] The parse position within \p str
            std::string::size_type& pos,
-            //! The fixed-format width to use (if comma is not found)
-            const std::string::size_type width,
-            //! The input is known to be free-format
-            const bool free_format = false
+            std::string::size_type len
        );


--- a/src/fileFormats/stl/STLCore.C
+++ b/src/fileFormats/stl/STLCore.C
@ -135,9 +135,8 @@ int Foam::fileFormats::STLCore::detectBinaryHeader

        bad =
            (
-                dataFileSize < STLHeaderSize
-             || nTris < (dataFileSize - STLHeaderSize)/50
-             || nTris > (dataFileSize - STLHeaderSize)/25
+                nTris < int(dataFileSize - STLHeaderSize)/50
+             || nTris > int(dataFileSize - STLHeaderSize)/25
            );
    }

@ -209,9 +208,8 @@ Foam::fileFormats::STLCore::readBinaryHeader

        bad =
            (
-                dataFileSize < STLHeaderSize
-             || nTris < (dataFileSize - STLHeaderSize)/50
-             || nTris > (dataFileSize - STLHeaderSize)/25
+                nTris < int(dataFileSize - STLHeaderSize)/50
+             || nTris > int(dataFileSize - STLHeaderSize)/25
            );
    }

--- a/src/finiteArea/Make/files
+++ b/src/finiteArea/Make/files
@ -55,6 +55,7 @@ $(basicFaPatchFields)/transform/transformFaPatchScalarField.C
 constraintFaPatchFields = $(faPatchFields)/constraint
 $(constraintFaPatchFields)/empty/emptyFaPatchFields.C
 $(constraintFaPatchFields)/processor/processorFaPatchFields.C
+$(constraintFaPatchFields)/processor/processorFaPatchScalarField.C
 $(constraintFaPatchFields)/wedge/wedgeFaPatchFields.C
 $(constraintFaPatchFields)/wedge/wedgeFaPatchScalarField.C
 $(constraintFaPatchFields)/cyclic/cyclicFaPatchFields.C
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -180,26 +180,12 @@ Foam::faBoundaryMesh::faBoundaryMesh

 void Foam::faBoundaryMesh::calcGeometry()
 {
-    // processor initGeometry send/recv the following:
-    //   - edgeCentres() : faMesh::edgeCentres()
-    //   - edgeLengths() : faMesh::Le()
-    //   - edgeFaceCentres() : faMesh::areaCentres()
-    //
-    // faMesh::Le() has its own point-to-point communication (OK) but
-    // triggers either/or edgeAreaNormals(), pointAreaNormals()
-    // with their own communication that can block.
-
+    // processorFaPatch geometry triggers calculation of pointNormals.
    // This uses parallel comms and hence will not be trigggered
    // on processors that do not have a processorFaPatch so instead
    // force construction.
-
-    (void)mesh_.edgeAreaNormals();
    (void)mesh_.pointAreaNormals();

-    (void)mesh_.areaCentres();
-    (void)mesh_.faceAreaNormals();
-
-
    PstreamBuffers pBufs(Pstream::defaultCommsType);

    if
@ -787,15 +773,12 @@ bool Foam::faBoundaryMesh::checkDefinition(const bool report) const

 void Foam::faBoundaryMesh::movePoints(const pointField& p)
 {
-    // See comments in calcGeometry()
-
-    (void)mesh_.edgeAreaNormals();
+    // processorFaPatch geometry triggers calculation of pointNormals.
+    // This uses parallel comms and hence will not be trigggered
+    // on processors that do not have a processorFaPatch so instead
+    // force construction.
    (void)mesh_.pointAreaNormals();

-    (void)mesh_.areaCentres();
-    (void)mesh_.faceAreaNormals();
-
-
    PstreamBuffers pBufs(Pstream::defaultCommsType);

    if
--- a/src/finiteArea/faMesh/faMesh.C
+++ b/src/finiteArea/faMesh/faMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2020-2023 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -277,27 +277,6 @@ void Foam::faMesh::clearOut() const

 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

-void Foam::faMesh::syncGeom()
-{
-    if (UPstream::parRun())
-    {
-        // areaCentres()
-        if (faceCentresPtr_)
-        {
-            faceCentresPtr_->boundaryFieldRef()
-                .evaluateCoupled<processorFaPatch>();
-        }
-
-        // faceAreaNormals()
-        if (faceAreaNormalsPtr_)
-        {
-            faceAreaNormalsPtr_->boundaryFieldRef()
-                .evaluateCoupled<processorFaPatch>();
-        }
-    }
-}
-
-
 bool Foam::faMesh::init(const bool doInit)
 {
    if (doInit)
@ -317,7 +296,18 @@ bool Foam::faMesh::init(const bool doInit)
    // Calculate the geometry for the patches (transformation tensors etc.)
    boundary_.calcGeometry();

-    syncGeom();
+    // Ensure processor/processor information is properly synchronised
+    if (Pstream::parRun())
+    {
+        const_cast<areaVectorField&>(areaCentres()).boundaryFieldRef()
+            .evaluateCoupled<processorFaPatch>();
+
+        // This roughly corresponds to what OpenFOAM-v2112 (and earlier) had,
+        // but should nominally be unnecessary.
+        //
+        /// const_cast<areaVectorField&>(faceAreaNormals()).boundaryFieldRef()
+        ///     .evaluateCoupled<processorFaPatch>();
+    }

    return false;
 }
@ -999,6 +989,7 @@ bool Foam::faMesh::movePoints()

    clearGeomNotAreas();

+    // To satisfy the motion interface for MeshObject, const cast is needed
    if (patchPtr_)
    {
        patchPtr_->movePoints(newPoints);
@ -1012,8 +1003,6 @@ bool Foam::faMesh::movePoints()

    // Note: Fluxes were dummy?

-    syncGeom();
-
    return true;
 }

--- a/src/finiteArea/faMesh/faMesh.H
+++ b/src/finiteArea/faMesh/faMesh.H
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -617,10 +617,6 @@ public:
        //- Initialise non-demand-driven data etc
        bool init(const bool doInit);

-        //- Processor/processor synchronisation for geometry fields.
-        //  Largely internal use only (slightly hacky).
-        void syncGeom();
-

    // Database

--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Kutalmis Bercin	55d5f99047	WIP: TUT: tutorial adjustment for backward compatibility	2023-01-31 16:14:43 +00:00
Kutalmis Bercin	d439440453	WIP: TUT: tutorial adjustments for backward compatibility	2023-01-20 16:56:11 +00:00
Andrew Heather	6c66cc58a6	WIP: Updated lagrangian solvers followinf library changes	2023-01-10 13:35:42 +00:00
Andrew Heather	265b3dd2f8	WIP: deprecated lagrangian coal combustion - models added to reactingMultiphase	2023-01-10 13:35:42 +00:00
Andrew Heather	b17a9128a1	ENH: Added new parcelFoam solver based on kinematicParcelFoam	2023-01-10 13:35:42 +00:00
Andrew Heather	31d4148173	STYLE: lagrangian turbulence library - removed 'basic' prefix	2023-01-10 13:35:42 +00:00
Andrew Heather	13b8a3c02a	STYLE: lagrangian intermediate library - removed 'basic' prefix	2023-01-10 13:35:42 +00:00
Andrew Heather	715831610d	WIP: Added new rhoParcelFoam solver based on reactingParcelFoam	2023-01-10 13:35:42 +00:00
Andrew Heather	ab42d1a1e6	WIP: migrated coalCloudList to parcelCloudModelList	2023-01-10 13:35:42 +00:00
Andrew Heather	13cf489f63	WIP: clouds - moved spray cloud and parcels to main intermediate library	2023-01-10 13:35:42 +00:00
Andrew Heather	3277fb931b	WIP: clouds - updated dependent code	2023-01-10 13:35:40 +00:00
Andrew Heather	9e83d31138	WIP: clouds - refactoring	2023-01-10 13:34:47 +00:00
Andrew Heather	878f1bf773	WIP: Clouds - new run-time selection based on openfoam.org changes	2023-01-10 13:34:47 +00:00
Andrew Heather	84348dc3f5	WIP: clouds - added debug - based on openfoam.org changes	2023-01-10 13:34:47 +00:00
Andrew Heather	1f7c038ad7	DUMMY	2023-01-10 13:34:47 +00:00
Andrew Heather	f20361d44a	WIP: Clouds - new run-time selection based on openfoam.org changes	2023-01-10 13:34:33 +00:00
Andrew Heather	580a003d4d	WIP: Removed base classes	2023-01-10 13:34:22 +00:00
mattijs	d74572cae8	CONFIG: set API level to 2212	2022-12-28 09:33:42 +00:00
Andrew Heather	1cd47839ac	CONFIG: Reset version post-release	2022-12-23 13:53:49 +00:00
				`@ -0,0 +1,2 @@`
				`auto parcels = parcelCloudModelList(g, rhoInf, U, muc);`
				`@ -0,0 +1 @@`
				`regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();`