WIP: duplicate constructor entries (#2538 )

WIP: Updated multiphaseEuler runtime constructor table names WIP: Updated multiphaseInter runtime constructor table names WIP: Updated reactingEuler runtime constructor table names WIP: Updated twoPhaseEuler runtime constructor table names
ENH: faceAgglomerate: read patch-based agglomeration. Fixes #2558 .
2022-08-19 11:01:42 +02:00 · 2022-08-11 11:12:55 +01:00 · 2022-08-09 15:58:40 +01:00 · 2022-08-08 12:50:23 +01:00 · 2022-08-04 17:09:39 +00:00 · 2022-08-04 17:09:38 +00:00
915 changed files with 8597 additions and 7248 deletions
--- a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
@ -116,7 +116,7 @@ int main(int argc, char *argv[])
            "normalisedGradP",
            tmagGradP()/max(tmagGradP())
        );
-        normalisedGradP.writeOpt() = IOobject::AUTO_WRITE;
+        normalisedGradP.writeOpt(IOobject::AUTO_WRITE);
        tmagGradP.clear();

        ++runTime;
--- a/applications/solvers/combustion/fireFoam/createClouds.H
+++ b/applications/solvers/combustion/fireFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-reactingCloud parcels
+basicReactingCloud 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 "reactingCloud.H"
+#include "basicReactingCloud.H"
 #include "surfaceFilmModel.H"
 #include "pyrolysisModelCollection.H"
 #include "radiationModel.H"
--- 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 "kinematicCloud.H"
-    #define kinematicTypeCloud kinematicCloud
+    #include "basicKinematicCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCloud
 #else
-    #include "kinematicCollidingCloud.H"
-    #define kinematicTypeCloud kinematicCollidingCloud
+    #include "basicKinematicCollidingCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCollidingCloud
 #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
-        kCloud.storeGlobalPositions();
+        kinematicCloud.storeGlobalPositions();

        mesh.update();

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

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

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

-        fvVectorMatrix cloudSU(kCloud.SU(Uc));
+        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
@ -12,8 +12,6 @@ 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,11 +11,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)/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 "kinematicCloud.H"
-    #define kinematicTypeCloud kinematicCloud
+    #include "basicKinematicCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCloud
 #else
-    #include "kinematicCollidingCloud.H"
-    #define kinematicTypeCloud kinematicCollidingCloud
+    #include "basicKinematicCollidingCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCollidingCloud
 #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 " << kCloud.name() << endl;
-        kCloud.evolve();
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+        kinematicCloud.evolve();

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

-        fvVectorMatrix cloudSU(kCloud.SU(Uc));
+        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
@ -10,8 +10,6 @@ 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,11 +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/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;
-kinematicTypeCloud kCloud
+basicKinematicTypeCloud kinematicCloud
 (
    kinematicCloudName,
-    g,
    rhoc,
    Uc,
-    muc
+    muc,
+    g
 );

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

 // Update alphac from the particle locations
-alphac = max(1.0 - kCloud.theta(), alphacMin);
+alphac = max(1.0 - kinematicCloud.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,6 +6,7 @@ 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 \
@ -34,6 +35,7 @@ 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 "thermoCloud.H"
-#include "reactingMultiphaseCloud.H"
+#include "basicThermoCloud.H"
+#include "coalCloud.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;
-reactingMultiphaseCloud coalParcels
+coalCloud coalParcels
 (
    "coalCloud1",
-    g,
    rho,
    U,
+    g,
    slgThermo
 );

 Info<< "\nConstructing limestone cloud" << endl;
-thermoCloud limestoneParcels
+basicThermoCloud limestoneParcels
 (
    "limestoneCloud1",
-    g,
    rho,
    U,
+    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
@ -9,8 +9,6 @@ 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;
-kinematicCollidingCloud kCloud
+basicKinematicCollidingCloud kinematicCloud
 (
    kinematicCloudName,
-    g,
    rhoInf,
    U,
-    mu
+    mu,
+    g
 );

 IOobject Hheader
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
@ -10,8 +10,6 @@ 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 "kinematicCollidingCloud.H"
+#include "basicKinematicCollidingCloud.H"

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

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

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

        mesh.update();

        U.correctBoundaryConditions();

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

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

-        kCloud.evolve();
+        kinematicCloud.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 "kinematicCollidingCloud.H"
+#include "basicKinematicCollidingCloud.H"

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

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

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

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

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

        runTime.write();

--- a/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/kinematicParcelFoam/Make/options
@ -8,12 +8,6 @@ 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)
     ==
-        scalar(1)/rhoInfValue*parcels.SU(U)
+        parcels.SU(U, true)
      + fvOptions(U)
    );

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

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

--- 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 "kinematicCloud.H"
+#include "basicKinematicCloud.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
 #include "CorrectPhi.H"
@ -67,6 +67,7 @@ int main(int argc, char *argv[])
    #include "createDyMControls.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
+    #include "createRegionControls.H"
    #include "createUfIfPresent.H"

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

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

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

-        if (pimple.solveFlow())
+        if (solvePrimaryRegion)
        {
            // --- PIMPLE loop
            while (pimple.loop())
--- a/applications/solvers/lagrangian/parcelFoam/Make/files
+++ b/applications/solvers/lagrangian/parcelFoam/Make/files
@ -1,3 +0,0 @@
-parcelFoam.C
-
-EXE = $(FOAM_USER_APPBIN)/parcelFoam
--- a/applications/solvers/lagrangian/parcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/parcelFoam/UEqn.H
@ -1,22 +0,0 @@
-    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
@ -1,2 +0,0 @@
-auto parcels = parcelCloudModelList(g, rhoInf, U, muc);
-
--- a/applications/solvers/lagrangian/parcelFoam/createFieldRefs.H
+++ b/applications/solvers/lagrangian/parcelFoam/createFieldRefs.H
@ -1 +0,0 @@
-regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
--- a/applications/solvers/lagrangian/parcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/parcelFoam/createFields.H
@ -1,85 +0,0 @@
-#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
@ -1,58 +0,0 @@
-
-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
@ -1,154 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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;
-reactingTypeCloud parcels
+basicReactingTypeCloud 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 reactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
+#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)

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

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

@ -104,7 +105,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 (pimple.solveFlow() && correctPhi)
+        if (solvePrimaryRegion && correctPhi)
        {
            divrhoU.reset
            (
@ -132,7 +133,7 @@ int main(int argc, char *argv[])

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

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

-        if (pimple.solveFlow())
+        if (solvePrimaryRegion)
        {
            if (pimple.nCorrPIMPLE() <= 1)
            {
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
@ -7,6 +7,7 @@ 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;
-reactingTypeCloud parcels
+basicReactingTypeCloud 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 reactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
+#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)

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

--- a/applications/solvers/lagrangian/rhoParcelFoam/EEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/EEqn.H
@ -1,41 +0,0 @@
-{
-    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
@ -1,3 +0,0 @@
-rhoParcelFoam.C
-
-EXE = $(FOAM_USER_APPBIN)/rhoParcelFoam
--- a/applications/solvers/lagrangian/rhoParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/rhoParcelFoam/UEqn.H
@ -1,34 +0,0 @@
-    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
@ -1,51 +0,0 @@
-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
@ -1,3 +0,0 @@
-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
@ -1,5 +0,0 @@
-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
@ -1,132 +0,0 @@
-#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
@ -1,35 +0,0 @@
-#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
@ -1,6 +0,0 @@
-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
@ -1,96 +0,0 @@
-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
@ -1,50 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
@ -1,206 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
@ -1,55 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
@ -1,137 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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,6 +7,7 @@ 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 \
@ -46,6 +47,7 @@ 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;
-reactingMultiphaseCloud parcels
+coalCloud 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 "reactingMultiphaseCloud.H"
+#include "coalCloud.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,6 +35,7 @@ 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;
-sprayCloud parcels
+basicSprayCloud parcels
 (
    "sprayCloud",
-    g,
    rho,
    U,
+    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
@ -39,6 +39,7 @@ 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 "sprayCloud.H"
+#include "basicSprayCloud.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,6 +34,7 @@ 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,6 +40,7 @@ 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 "sprayCloud.H"
+#include "basicSprayCloud.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 "sprayCloud.H"
+#include "basicSprayCloud.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,8 +9,6 @@ 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;
-kinematicCloud kCloud
+basicKinematicCloud kinematicCloud
 (
    kinematicCloudName,
-    g,
    rho,
    U,
-    thermo.mu()
+    thermo.mu(),
+    g
 );
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
@ -10,8 +10,6 @@ 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 "kinematicCloud.H"
+#include "basicKinematicCloud.H"

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

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

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

        mesh.update();

        U.correctBoundaryConditions();

-        Info<< "Evolving " << kCloud.name() << endl;
-        kCloud.evolve();
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+        kinematicCloud.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 "kinematicCloud.H"
+#include "basicKinematicCloud.H"

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

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

-        Info<< "Evolving " << kCloud.name() << endl;
-        kCloud.evolve();
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+        kinematicCloud.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 "kinematicCloud.H"
+#include "basicKinematicCloud.H"

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

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

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

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

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

        // Update continuous phase volume fraction field
-        alphac = max(1.0 - kCloud.theta(), alphacMin);
+        alphac = max(1.0 - kinematicCloud.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(kCloud.SU(U));
+        fvVectorMatrix cloudSU(kinematicCloud.SU(U));
        volVectorField cloudVolSUSu
        (
            IOobject
--- a/applications/solvers/multiphase/MPPICInterFoam/Make/options
+++ b/applications/solvers/multiphase/MPPICInterFoam/Make/options
@ -12,14 +12,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/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;
-kinematicCloud kCloud
+basicKinematicCloud kinematicCloud
 (
    "kinematicCloud",
-    g,
    rho,
    U,
-    mu
+    mu,
+    g
 );

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

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

 surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleIO.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleIO.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -59,7 +59,7 @@ Foam::DTRMParticle::DTRMParticle
 {
    if (readFields)
    {
-        if (is.format() == IOstream::ASCII)
+        if (is.format() == IOstreamOption::ASCII)
        {
            is >> p0_ >> p1_ >> I0_ >> I_ >> dA_ >> transmissiveId_;
        }
@ -115,7 +115,7 @@ void Foam::DTRMParticle::writeProperties

 Foam::Ostream& Foam::operator<<(Ostream& os, const DTRMParticle& p)
 {
-    if (os.format() == IOstream::ASCII)
+    if (os.format() == IOstreamOption::ASCII)
    {
        os  << static_cast<const particle&>(p)
            << token::SPACE << p.p0_
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C
@ -699,7 +699,7 @@ void Foam::radiation::laserDTRM::calculate()

            for (label pointi = 0; pointi < lines.size(); pointi += 2)
            {
-                os.write(linePointRef(lines[pointi], lines[pointi+1]));
+                os.writeLine(lines[pointi], lines[pointi+1]);
            }
        }

--- a/applications/test/CompactIOList/Test-CompactIOList.C
+++ b/applications/test/CompactIOList/Test-CompactIOList.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -49,8 +49,8 @@ int main(int argc, char *argv[])
    #include "setRootCase.H"
    #include "createTime.H"

-    IOstream::streamFormat format = IOstream::BINARY;
-    // IOstream::streamFormat format = IOstream::ASCII;
+    IOstreamOption streamOpt(IOstreamOption::BINARY);
+    // IOstreamOption streamOpt(IOstreamOption::ASCII);

    const label size = 20000000;

@ -85,11 +85,7 @@ int main(int argc, char *argv[])
            << runTime.cpuTimeIncrement() << " s" << nl << endl;


-        faces2.writeObject
-        (
-            IOstreamOption(format),
-            true
-        );
+        faces2.writeObject(streamOpt, true);

        Info<< "Written old format faceList in = "
            << runTime.cpuTimeIncrement() << " s" << nl << endl;
@ -145,11 +141,7 @@ int main(int argc, char *argv[])
            << runTime.cpuTimeIncrement() << " s" << nl << endl;


-        faces2.writeObject
-        (
-            IOstreamOption(format),
-            true
-        );
+        faces2.writeObject(streamOpt, true);

        Info<< "Written new format faceList in = "
            << runTime.cpuTimeIncrement() << " s" << nl << endl;
--- a/applications/test/IndirectList/Test-IndirectList.C
+++ b/applications/test/IndirectList/Test-IndirectList.C
@ -146,7 +146,7 @@ int main(int argc, char *argv[])
    {
        Info<<"Writing output to " << binaryOutput << endl;

-        OFstream os(binaryOutput, IOstream::BINARY);
+        OFstream os(binaryOutput, IOstreamOption::BINARY);

        os.writeEntry("idl1", idl1);
        os.writeEntry("idl2", idl2);
--- a/applications/test/TestTools/TestTools.H
+++ b/applications/test/TestTools/TestTools.H
@ -28,15 +28,14 @@ Description

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

-using namespace Foam;
-
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "Matrix.H"
 #include "Random.H"
 #include <chrono>

+using namespace Foam;
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 // Total number of unit tests
--- a/applications/test/decomposedBlockData/Test-decomposedBlockData.C
+++ b/applications/test/decomposedBlockData/Test-decomposedBlockData.C
@ -78,8 +78,8 @@ int main(int argc, char *argv[])
    OFstream os
    (
        objPath,
-        IOstream::BINARY,
-        IOstream::currentVersion,
+        IOstreamOption::BINARY,
+        IOstreamOption::currentVersion,
        runTime.writeCompression()
    );
    if (!os.good())
--- a/applications/test/dimField/Test-dimField.C
+++ b/applications/test/dimField/Test-dimField.C
@ -42,7 +42,7 @@ namespace Foam
    // Something like an internal state field. Probably only dimensionless
    typedef DimensionedField<uint8_t, volMesh> dimUint8Field;

-    defineTemplate2TypeNameAndDebug(dimUint8Field, 0);
+    defineTemplateTypeNameAndDebug(dimUint8Field, 0);

 } // End namespace Foam
 #endif
--- a/applications/test/edges/Test-edges.C
+++ b/applications/test/edges/Test-edges.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -82,8 +82,7 @@ int main(int argc, char *argv[])

    Info<< nl << "hash-like functionality" << nl;

-    // doesn't work e4 = -1;
-    e4.start() = e4.end() = -1;
+    e4.clear();

    printInfo(e4);
    for (label i : {2, -1, 2, 1, 4, 1, 2, 3})
@ -93,24 +92,24 @@ int main(int argc, char *argv[])
        printInfo(e4);
    }

-    e4.start() = e4.end() = -1;
+    e4.clear();
    Info<< "insert from list\n";
    labelHashSet newIndices({2, -1, 2, 1, 4, 1, 2, 3});
    e4.insert(newIndices.toc());
    printInfo(e4);

-    e4.start() = e4.end() = -1;
+    e4.clear();
    Info<< "insert from list\n";
    e4.insert({0, 5, 2, -1, 2, 1, 4, 1, 2, 3});
    printInfo(e4);

    FixedList<label, 8> otherIndices{12, 2, -1, 1, 4, 1, 2, 3};
-    e4.start() = e4.end() = -1;
+    e4.clear();
    Info<< "insert from list: " << otherIndices << nl;
    e4.insert(otherIndices);
    printInfo(e4);

-    e4.start() = e4.end();
+    e4.a() = e4.b();
    Info<< "erase from list: " << otherIndices << nl;
    Info<< "removed " << e4.erase(otherIndices) << " values" << nl;
    printInfo(e4);
--- a/applications/test/foamToMetisGraph/foamToMetisGraph.C
+++ b/applications/test/foamToMetisGraph/foamToMetisGraph.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -56,7 +56,7 @@ int main(int argc, char *argv[])
    // No. of Nodes = nCells
    // No. of Edges connecting Nodes = nInternalFaces

-    OFstream os(args.caseName() + ".graph", IOstream::ASCII);
+    OFstream os(args.caseName() + ".graph", IOstreamOption::ASCII);

    os  << "%% metis graph file, of an OpenFOAM mesh %%" << nl
        << "%% nCells=" << mesh.nCells()
--- a/applications/test/fvc2D/Test-fvc2D.C
+++ b/applications/test/fvc2D/Test-fvc2D.C
@ -39,7 +39,7 @@ namespace Foam
 {
    typedef GeometricField<vector2D, fvPatchField, volMesh> volVector2DField;

-    defineTemplate2TypeNameAndDebug(volVector2DField::Internal, 0);
+    defineTemplateTypeNameAndDebug(volVector2DField::Internal, 0);
    defineTemplateTypeNameAndDebug(volVector2DField, 0);

    typedef fvPatchField<vector2D> fvPatchVector2DField;
--- a/applications/test/instant/Test-instant.C
+++ b/applications/test/instant/Test-instant.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -30,11 +30,31 @@ Description

 #include "argList.H"
 #include "instant.H"
+#include "Pair.H"
 #include "fileNameInstant.H"
 #include "DynamicList.H"

 using namespace Foam;

+template<class T>
+Ostream& printInstant(const UList<T>& times, const label i)
+{
+    if (i >= 0 && i < times.size())
+    {
+        Info<< " (" << times[i] << ")";
+    }
+    return Info;
+}
+
+template<class T>
+Ostream& printInstant(const UList<T>& times, const Pair<label>& range)
+{
+    printInstant(times, range.first());
+    printInstant(times, range.second());
+    return Info;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:

@ -48,6 +68,7 @@ int main(int argc, char *argv[])
    times.append({300.456, "def"});
    times.append({454.456, "xyz"});
    times.append({10, "ten"});
+    times.append({15, "fifteen"});

    {
        word timeName("twenty");
@ -61,6 +82,19 @@ int main(int argc, char *argv[])
    sort(times);
    Info<< "Sorted:" << times << nl;

+    for (const scalar val : { -0.5, 5.0, 18.0, 25.0, 450.0, 480.0 })
+    {
+        label start = instant::findStart(times, val);
+        Pair<label> range = instant::findRange(times, val);
+
+        Info<< nl
+            << "time:" << val << nl;
+        Info<< "    start:" << start;
+        printInstant(times, start) << nl;
+
+        Info<< "    range:" << range;
+        printInstant(times, range) << nl;
+    }

    DynamicList<fileNameInstant> files;
    files.append(fileNameInstant{});
--- a/applications/test/limits/Test-limits.C
+++ b/applications/test/limits/Test-limits.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -24,10 +24,11 @@ License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

 Description
-    Print max limits.
+    Print some numerical limits.

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

+#include <cmath>
 #include <limits>
 #include "int.H"
 #include "uint.H"
@ -37,21 +38,63 @@ Description

 using namespace Foam;

+std::ostream& print(const char* tag, float val)
+{
+    std::cout
+        << tag << val
+        << " 0x" << std::hex << *reinterpret_cast<const uint32_t*>(&val);
+    return std::cout;
+}
+
+
+std::ostream& print(const char* tag, double val)
+{
+    std::cout
+        << tag << val
+        << " 0x" << std::hex << *reinterpret_cast<const uint64_t*>(&val);
+    return std::cout;
+}
+
+
+// Have (float|double)Scalar(GREAT|SMALL|..)
+
+#define PrintFloatLimits(FloatType, NanFunction)                              \
+{                                                                             \
+    print("max:", std::numeric_limits<FloatType>::max());                     \
+    print(" VGREAT:", FloatType##ScalarVGREAT);                               \
+    print(" ROOTVGREAT:", FloatType##ScalarROOTVGREAT) << nl;                 \
+                                                                              \
+    print("min:", std::numeric_limits<FloatType>::min());                     \
+    print(" VSMALL:", FloatType##ScalarVSMALL);                               \
+    print(" ROOTVSMALL:", FloatType##ScalarROOTVSMALL) << nl;                 \
+                                                                              \
+    print("epsilon:", std::numeric_limits<FloatType>::epsilon());             \
+    print(" SMALL:", FloatType##ScalarSMALL);                                 \
+    print(" ROOTSMALL:", FloatType##ScalarROOTSMALL) << nl;                   \
+                                                                              \
+    print("1/epsilon:", 1/std::numeric_limits<FloatType>::epsilon());         \
+    print(" GREAT:", FloatType##ScalarGREAT);                                 \
+    print(" ROOTGREAT:", FloatType##ScalarROOTGREAT) << nl;                   \
+                                                                              \
+    print("nan:", std::NanFunction("nan")) << nl;                             \
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:

 int main(int argc, char *argv[])
 {
-    //NONE Info<<"int16:" << pTraits<int16_t>::max << nl;
-    Info<< "=max=" << nl;
-    Info<< "uint8:" << pTraits<uint8_t>::max << nl;
-    Info<< "int16:" << std::numeric_limits<int16_t>::max() << nl;
-    Info<< "int32:" << pTraits<int32_t>::max << nl;
-    Info<< "uint32:" << pTraits<uint32_t>::max << nl;
-    Info<< "int64:" << pTraits<int64_t>::max << nl;
-    Info<< "uint64:" << pTraits<uint64_t>::max << nl;
+    //NONE cout<<"int16:" << pTraits<int16_t>::max << nl;
+    cout<< "=max=" << nl;
+    cout<< "uint8:" << pTraits<uint8_t>::max << nl;
+    cout<< "int16:" << std::numeric_limits<int16_t>::max() << nl;
+    cout<< "int32:" << pTraits<int32_t>::max << nl;
+    cout<< "uint32:" << pTraits<uint32_t>::max << nl;
+    cout<< "int64:" << pTraits<int64_t>::max << nl;
+    cout<< "uint64:" << pTraits<uint64_t>::max << nl;

-    Info<< nl;
+    cout<< nl;

    cout<< "int16:" << std::numeric_limits<int16_t>::max() << nl;
    cout<< "int32:" << pTraits<int32_t>::max << nl;
@ -59,7 +102,7 @@ int main(int argc, char *argv[])
    cout<< "int64:" << pTraits<int64_t>::max << nl;
    cout<< "uint64:" << pTraits<uint64_t>::max << nl;

-    Info<< nl << "=digits=" << nl;
+    cout<< nl << "=digits=" << nl;

    cout<< "int16:" << std::numeric_limits<int16_t>::digits << nl;
    cout<< "int32:" << std::numeric_limits<int32_t>::digits << nl;
@ -70,29 +113,15 @@ int main(int argc, char *argv[])
    cout<< "double:" << std::numeric_limits<double>::digits << nl;
    cout<< "long double:" << std::numeric_limits<long double>::digits << nl;

-    Info<< nl << "=float=" << nl;
+    // std::nanl (long double)
+    cout<< nl;
+    cout<< nl << "=float=" << nl;
+    PrintFloatLimits(float, nanf);

-    cout<< "max:" << std::numeric_limits<float>::max()
-        << " VGREAT:" << floatScalarVGREAT << nl;
-    cout<< "min:" << std::numeric_limits<float>::min()
-        << " VSMALL:" << floatScalarVSMALL << nl;
-    cout<< "epsilon:" << std::numeric_limits<float>::epsilon()
-        << " SMALL:" << floatScalarSMALL << nl;
-    cout<< "1/epsilon:" << 1.0f/std::numeric_limits<float>::epsilon()
-        << " GREAT:" << floatScalarGREAT << nl;
+    cout<< nl << "=double=" << nl;
+    PrintFloatLimits(double, nan);

-    Info<< nl << "=double=" << nl;
-
-    cout<< "max:" << std::numeric_limits<double>::max()
-        << " VGREAT:" << doubleScalarVGREAT << nl;
-    cout<< "min:" << std::numeric_limits<double>::min()
-        << " VSMALL:" << doubleScalarVSMALL << nl;
-    cout<< "epsilon:" << std::numeric_limits<double>::epsilon()
-        << " SMALL:" << doubleScalarSMALL << nl;
-    cout<< "1/epsilon:" << 1.0f/std::numeric_limits<double>::epsilon()
-        << " GREAT:" << doubleScalarGREAT << nl;
-
-    Info << "---\nEnd\n" << endl;
+    cout<< "---\nEnd\n" << std::endl;

    return 0;
 }
--- a/applications/test/matrices/DiagonalMatrix/Test-DiagonalMatrix.C
+++ b/applications/test/matrices/DiagonalMatrix/Test-DiagonalMatrix.C
@ -38,8 +38,7 @@ Description

 #include "DiagonalMatrix.H"
 #include "RectangularMatrix.H"
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "TestTools.H"

--- a/applications/test/matrices/RectangularMatrix/Test-RectangularMatrix.C
+++ b/applications/test/matrices/RectangularMatrix/Test-RectangularMatrix.C
@ -44,8 +44,7 @@ Note

 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "TestTools.H"
--- a/applications/test/matrices/SquareMatrix/Test-SquareMatrix.C
+++ b/applications/test/matrices/SquareMatrix/Test-SquareMatrix.C
@ -45,8 +45,7 @@ Note
 #include "scalarMatrices.H"
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "TestTools.H"
--- a/applications/test/matrices/SymmetricSquareMatrix/Test-SymmetricSquareMatrix.C
+++ b/applications/test/matrices/SymmetricSquareMatrix/Test-SymmetricSquareMatrix.C
@ -46,8 +46,7 @@ Note
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
 #include "SymmetricSquareMatrix.H"
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "Random.H"
--- a/applications/test/primitives/Test-primitives.C
+++ b/applications/test/primitives/Test-primitives.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -334,7 +334,7 @@ int main(int argc, char *argv[])

        DynamicList<char> buf;

-        OListStream os(std::move(buf), IOstream::BINARY);
+        OListStream os(std::move(buf), IOstreamOption::BINARY);
        os << srcList;

        os.swap(buf); // Recover buffer
@ -342,7 +342,7 @@ int main(int argc, char *argv[])
        // Read back
        List<scalar> dstList;

-        UIListStream is(buf, IOstream::BINARY);
+        UIListStream is(buf, IOstreamOption::BINARY);
        is.setScalarByteSize(sizeof(otherType));

        Info<< "Stream scalar-size ("
--- a/applications/test/volField/Test-volField.C
+++ b/applications/test/volField/Test-volField.C
@ -42,7 +42,7 @@ namespace Foam
    // - still needs some basic boundary conditions!!
    typedef GeometricField<uint8_t, fvPatchField, volMesh> volUint8Field;

-    defineTemplate2TypeNameAndDebug(volUint8Field, 0);
+    defineTemplateTypeNameAndDebug(volUint8Field, 0);

 } // End namespace Foam
 #endif
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
@ -1397,26 +1397,20 @@ void extrudeGeometricProperties
        {
            Pout<< "Model     :" << faceCentres[facei] << endl
                << "regionMesh:" << regionMesh.faceCentres()[facei] << endl;
-            faceStr.write
+            faceStr.writeLine
            (
-                linePointRef
-                (
-                    faceCentres[facei],
-                    regionMesh.faceCentres()[facei]
-                )
+                faceCentres[facei],
+                regionMesh.faceCentres()[facei]
            );
        }
        forAll(cellCentres, celli)
        {
            Pout<< "Model     :" << cellCentres[celli] << endl
                << "regionMesh:" << regionMesh.cellCentres()[celli] << endl;
-            cellStr.write
+            cellStr.writeLine
            (
-                linePointRef
-                (
-                    cellCentres[celli],
-                    regionMesh.cellCentres()[celli]
-                )
+                cellCentres[celli],
+                regionMesh.cellCentres()[celli]
            );
        }
    }
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/PrintTable/PrintTable.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/PrintTable/PrintTable.C
@ -152,7 +152,7 @@ void Foam::PrintTable<KeyType, DataType>::print

        os  << endl;

-        const List<KeyType>& sortedTable = combinedTable.sortedToc();
+        const List<KeyType> sortedTable(combinedTable.sortedToc());

        forAll(sortedTable, keyI)
        {
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@ -1687,7 +1687,7 @@ void Foam::conformalVoronoiMesh::move()
                    )
                )
                {
-                    multipleIntersections.write(linePointRef(ptA, ptB));
+                    multipleIntersections.writeLine(ptA, ptB);
                }
            }
        }
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@ -83,7 +83,7 @@ SourceFiles
 namespace Foam
 {

-// Forward declaration of classes
+// Forward Declarations
 class initialPointsMethod;
 class relaxationModel;
 class faceAreaWeightModel;
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@ -1433,18 +1433,15 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 //                         && (mag(snapDir & norm[0]) > 0.5)
 //                        )
 //                        {
-//                            snapping1.write
+//                            snapping1.writeLine
 //                            (
-//                                linePointRef(dual, nearestPointOnTet)
+//                                dual,
+//                                nearestPointOnTet
 //                            );
-//
-//                            snapping2.write
+//                            snapping2.writeLine
 //                            (
-//                                linePointRef
-//                                (
-//                                    nearestPointOnTet,
-//                                    hitInfo.hitPoint()
-//                                )
+//                                nearestPointOnTet,
+//                                hitInfo.hitPoint()
 //                            );
 //
 //                            pts[cit->cellIndex()] = hitInfo.hitPoint();
@ -1764,9 +1761,10 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches

                    if ((*vcit)->real())
                    {
-                        featurePointDualsStr.write
+                        featurePointDualsStr.writeLine
                        (
-                            linePointRef(topoint(vit->point()), (*vcit)->dual())
+                            topoint(vit->point()),
+                            (*vcit)->dual()
                        );
                    }
                }
@ -1867,7 +1865,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                        << " " << vc2->dual()
                        << endl;

-                    startCellStr.write(linePointRef(vc1->dual(), vc2->dual()));
+                    startCellStr.writeLine(vc1->dual(), vc2->dual());

                    // Get patch by getting face between cells and the two
                    // points on the face that are not the feature vertex
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@ -1395,10 +1395,7 @@ void Foam::conformalVoronoiMesh::writePointPairs

        if (ptPairs_.isPointPair(vA, vB))
        {
-            os.write
-            (
-                linePointRef(topoint(vA->point()), topoint(vB->point()))
-            );
+            os.writeLine(topoint(vA->point()), topoint(vB->point()));
        }
    }
 }
--- a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -73,7 +73,7 @@ scalar getMergeDistance
    Info<< "Merge tolerance : " << mergeTol << nl
        << "Write tolerance : " << writeTol << endl;

-    if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
+    if (runTime.writeFormat() == IOstreamOption::ASCII && mergeTol < writeTol)
    {
        FatalErrorInFunction
            << "Your current settings specify ASCII writing with "
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@ -616,7 +616,7 @@ scalar getMergeDistance
        << endl;

    // check writing tolerance
-    if (mesh.time().writeFormat() == IOstream::ASCII && !dryRun)
+    if (mesh.time().writeFormat() == IOstreamOption::ASCII && !dryRun)
    {
        const scalar writeTol = std::pow
        (
@ -1342,6 +1342,17 @@ int main(int argc, char *argv[])
    const snapParameters snapParams(snapDict, dryRun);


+    Info<< "Setting refinement level of surface to be consistent"
+        << " with curvature." << endl;
+    surfaces.setCurvatureMinLevelFields
+    (
+        refineParams.curvature(),
+        meshRefiner.meshCutter().level0EdgeLength()
+    );
+    Info<< "Checked curvature refinement in = "
+        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+
+

    // Add all the cellZones and faceZones
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--- a/applications/utilities/miscellaneous/OSspecific/addr2line/addr2line.C
+++ b/applications/utilities/miscellaneous/OSspecific/addr2line/addr2line.C
@ -43,7 +43,7 @@ Description
 static void usage();
 static void version();
 static std::string getLine(const std::string&, const std::string&);
-static std::string pOpen(const std::string&, int line=0);
+static std::string pipeOpen(const std::string& cmd, const int lineNum = 0);


 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -132,51 +132,50 @@ void version()
 }


-std::string pOpen(const std::string& cmd, int line)
+// Read up to and including lineNum from the piped command
+// Return the final line read
+std::string pipeOpen(const std::string& cmd, int lineNum)
 {
-    std::string res;
+    std::string str;

-    FILE* cmdPipe = popen(cmd.c_str(), "r");
-    if (!cmdPipe) return res;
+    FILE* handle = popen(cmd.c_str(), "r");
+    if (!handle) return str;

    char* buf = nullptr;
+    size_t len = 0;
+    ssize_t nread;

-    // Read line number of lines
-    for (int cnt = 0; cnt <= line; ++cnt)
+    // Read lineNum number of lines
+    for
+    (
+        int cnt = 0;
+        cnt <= lineNum && (nread = ::getline(&buf, &len, handle)) >= 0;
+        ++cnt
+    )
    {
-        size_t linecap = 0;
-        ssize_t linelen = ::getline(&buf, &linecap, cmdPipe);
-
-        if (linelen < 0)
+        if (cnt == lineNum)
        {
-            break;
-        }
+            // Retain the last line, trimming trailing newline
+            str.assign(buf);

-        if (cnt == line)
-        {
-            res = std::string(buf);
-
-            // Trim trailing newline
-            if (res.size())
+            if (str.size())
            {
-                res.resize(res.size()-1);
+                str.resize(str.size()-1);
            }
-            break;
        }
    }

-    if (buf) free(buf);
+    free(buf);
+    pclose(handle);

-    pclose(cmdPipe);
-
-    return res;
+    return str;
 }


 std::string getLine(const std::string& filename, const std::string& addr)
 {
    std::string line =
-        pOpen
+        pipeOpen
        (
            "echo 'image lookup -va " + addr
          + "'"
--- a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C
+++ b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -160,7 +160,7 @@ bool writeZones
        // Force writing as ASCII
        writeOk = meshObject.regIOobject::writeObject
        (
-            IOstreamOption(IOstream::ASCII, compression),
+            IOstreamOption(IOstreamOption::ASCII, compression),
            true
        );
    }
@ -355,7 +355,7 @@ int main(int argc, char *argv[])
            runTime
        );

-        if (runTime.writeFormat() == IOstream::ASCII)
+        if (runTime.writeFormat() == IOstreamOption::ASCII)
        {
            // Only do zones when converting from binary to ascii
            // The other way gives problems since working on dictionary level.
--- a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
+++ b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
@ -663,7 +663,11 @@ int main(int argc, char *argv[])
        Info<< "Merge tolerance : " << mergeTol << nl
            << "Write tolerance : " << writeTol << endl;

-        if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
+        if
+        (
+            runTime.writeFormat() == IOstreamOption::ASCII
+         && mergeTol < writeTol
+        )
        {
            FatalErrorInFunction
                << "Your current settings specify ASCII writing with "
--- a/applications/utilities/parallelProcessing/redistributePar/passivePositionParticle.H
+++ b/applications/utilities/parallelProcessing/redistributePar/passivePositionParticle.H
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -159,7 +159,7 @@ public:
            p.origProc = ppi.origProc();
            p.origId = ppi.origId();

-            if (os.format() == IOstream::ASCII)
+            if (os.format() == IOstreamOption::ASCII)
            {
                os  << p.position
                    << token::SPACE << p.celli
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
@ -232,13 +232,13 @@ vtk::outputOptions getOutputOptions(const argList& args)

        if (!args.found("ascii"))
        {
-            if (sizeof(floatScalar) != 4 || sizeof(label) != 4)
+            if (sizeof(float) != 4 || sizeof(label) != 4)
            {
                opts.ascii(true);

                WarningInFunction
                    << "Using ASCII rather than legacy binary VTK format since "
-                    << "floatScalar and/or label are not 4 bytes in size."
+                    << "float and/or label are not 4 bytes in size."
                    << nl << endl;
            }
            else
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Andrew Heather	742d1dc604	WIP: duplicate constructor entries (#2538 ) WIP: Updated multiphaseEuler runtime constructor table names WIP: Updated multiphaseInter runtime constructor table names WIP: Updated reactingEuler runtime constructor table names WIP: Updated twoPhaseEuler runtime constructor table names	2022-08-19 11:01:42 +02:00
mattijs	1e02a4ae92	ENH: faceAgglomerate: read patch-based agglomeration. Fixes #2558 . Read from optional subdictionary.	2022-08-11 11:12:55 +01:00
Kutalmis Bercin	f249022bfe	ENH: binModel: make decomposePatchValues entry optional (#2530 ) DOC: forceCoeffs: correct the header file content (#2530)	2022-08-09 15:58:40 +01:00
mattijs	c418c28c66	ENH: snappyHexMesh : refine based on curvature See https://develop.openfoam.com/Development/openfoam/-/merge_requests/558 Fix up of per-region specificiation.	2022-08-08 12:50:23 +01:00
Mattijs Janssens	5ea365a2be	Merge branch 'feature-ep_1823_curvature_refinement' into 'develop' snappyHexMesh : refine based on curvature See merge request Development/openfoam!558	2022-08-04 17:09:39 +00:00
Mattijs Janssens	27c3d0c23b	snappyHexMesh : refine based on curvature	2022-08-04 17:09:38 +00:00
Mattijs Janssens	227727d413	Merge branch 'feature-vf-ext' into 'develop' Update of view factor generation using 2AI and 2LI methods plus CGAL for ray tracing See merge request Development/openfoam!551	2022-08-04 14:18:23 +00:00
mattijs	c08fc5ecd9	ENH: viewFactorGen: revert to v2206 without CGAL	2022-08-04 14:59:12 +01:00
sergio	994addc543	ENH: Changing key entry name to GaussQuadTol	2022-08-04 14:59:12 +01:00
mattijs	f2f71f6847	ENH: viewFactor: compile without lib	2022-08-04 14:59:12 +01:00
sergio	c652af4b82	STY: Minor style changes	2022-08-04 14:59:12 +01:00
sergio	bfef08a89f	ENH: Using globalIndex to create full triSurface	2022-08-04 14:59:12 +01:00
sergio	3b40ee8f6b	ENH: Updating tutorials for view factors generation	2022-08-04 14:59:12 +01:00
sergio	f5598706c5	ENH: Making agglomeration optional for viewFactor model	2022-08-04 14:59:11 +01:00
sergio	457979a7b7	ENH: Adding new viewFactor utility using CGAL	2022-08-04 14:59:11 +01:00
sergio	565c68f454	ENH: Adding intersection margen	2022-08-04 14:59:11 +01:00
sergio	b9507c21f9	ENH: Adding viewFactorsGenExt which uses pbrt for ray tracing	2022-08-04 14:59:11 +01:00
Takuya Yamamoto	88da2d5877	ENH: interfaceProperties: add smoother for interfacial curvatures (closes #2531 ) The new algorithm introduces an optional entry 'nAlphaSmoothCurvature', and aims to smooth interface curvatures to reduce spurious currents.	2022-08-04 14:39:33 +01:00
mattijs	77ecc4c6e0	BUG: sortedToc: reference to copy. Fixes #2554	2022-08-03 15:05:25 +01:00
Kutalmis Bercin	542fc0c1d1	ENH: lnGradSchemes: new uncorrectLnGrad scheme	2022-08-03 13:28:51 +01:00
Kutalmis Bercin	9ccc2d8fd5	ENH: edgeInterpolation: avoid division-by-zero errors in skew corrections	2022-08-03 13:28:44 +01:00
mattijs	c2cae92fc5	ENH: changeDictionary: support collated format. Fixes #2533	2022-08-01 15:16:28 +01:00
Kutalmis Bercin	5894874920	TUT: faSchemes: replace snGradSchemes with lnGradSchemes (fixes #2543 ) When a finite-area case could not find an entry for "lnGradSchemes" in the "faSchemes" file, the "corrected" scheme has been picked up by default. Therefore, any changes in "snGradSchemes" entry will not be read by finite-area models.	2022-07-25 15:06:51 +01:00
Mark Olesen	6320bab2b5	STYLE: IOstreams with float/double instead of floatScalar/doubleScalar - consistent with defining IO of int32_t/int64_t and with recent changes to ensightFile. Using the primitives directly instead of typedefs to them makes the code somewhat less opaque.	2022-07-22 15:43:37 +02:00
Mark Olesen	a4a8f77b7b	STYLE: more consistent use of CGAL_LIBS for surfaceBooleanFeatures	2022-07-21 17:23:07 +02:00
Mark Olesen	64f4745277	STYLE: various 'defineTemplate2...' macros are redundant - they are/were identical to versions without an embedded '2'	2022-07-21 17:07:16 +02:00
Mark Olesen	ac83b41aaf	ENH: improve demangled symbol names for safePrintStack - parse out symbols and use abi::__cxa_demangle for more readable names in safePrintStack. - shorten prefixed /path/openfoam/platforms/lib/... to start with "platforms/lib/..." to avoid unreadably long lines. - improved file-scope localization of helper functions. STYLE: use std::ios_base::basefield instead of dec\|oct\|hex for masking	2022-07-21 11:29:49 +02:00
mattijs	20f1afd9f7	ENH: tutorial: add comment	2022-07-20 13:38:31 +01:00
mattijs	873f7aac2d	BUG: multi-world: avoid on-the-fly comms. Fixes #2529 Probably can get away with less tetBasePtIs() triggering	2022-07-20 13:38:31 +01:00
mattijs	66ddf0a104	BUG: multi-world: only switch on warnComm if needed. See #2529	2022-07-20 13:38:30 +01:00
Mark Olesen	42150cf3fe	ENH: support output of face area-normals for boundaryData writer	2022-07-19 11:18:44 +02:00
Mark Olesen	5630db5493	ENH: robuster handling of mapMethod naming (#2535 ) - align timeVaryingMappedFixedValuePointPatchField keywords with MappedFile STYLE: minor cleanup of pointToPointPlanarInterpolation BUG: incorrect keyword for timeVaryingMappedFixedValuePointPatchField - lookup should be "fieldTable" (not "fieldTableName") for consistency with the output and other BCs. (Bug introduced by `a623ab42a3`)	2022-07-19 11:17:52 +02:00
Mark Olesen	dfdbe7efd0	ENH: add instance searching routines - find start index - find index range spanning a time	2022-07-19 11:17:52 +02:00
Mark Olesen	b4612b4c04	ENH: cleanup ensight surface reader (#2535 ) - some central (core) bits under fileFormats, - general surface reading relocated from sampling to surfMesh since it does not use any sampling-specific components and will permit re-use in meshTools (for example) - remove old mask, subDir methods from ensightFile which were previously relocated to ensightCase - improve handling of 'undef' values when generating and reading, respect Ensight component ordering when reading.	2022-07-19 11:17:52 +02:00
Mark Olesen	c4d18e97a3	ENH: additional methods for OBJstream - write point fields - writeLine (takes two points) - writeFace (takes list of face loop points)	2022-07-19 11:17:51 +02:00
Mark Olesen	dea31e9b4a	ENH: consistent member access for triFace / triangle etc. - can access the vertices/points as a(), b(), c()	2022-07-19 11:17:51 +02:00
Mark Olesen	3d892ace29	STYLE: set readOpt(..), writeOpt(..) by parameter, not by assignment STYLE: qualify format/version/compression with IOstreamOption not IOstream STYLE: reduce number of lookups when scanning {fa,fv}Solution STYLE: call IOobject::writeEndDivider as static	2022-07-19 11:17:47 +02:00
Mark Olesen	8b1514c99b	GIT: relocate rawIOField to src/OpenFOAM (closer to IOField classes) STYLE: RunFunctions - missed message for failure to find nFaces	2022-07-13 19:23:27 +02:00
Mark Olesen	1afd27db6c	BUG: incorrect ensight filename references (fixes #2532 ) - Ensight places restrictions both on variable names and on file names. When generating the variable to file name correspondence for use in the Ensight case file, previously used the less stringent variable name for both sides of the variable table. This would lead to situations where the (valid) variable name referred to the wrong file name. Now apply the file-name restriction consistently when creating the variable table. This is especially necessary since the stem of the filename additionally has specific characters (eg, ":<>[]") that can be problematic for the shell or file-system. ENH: avoid repeated '_' in qualified ensight names. - when replacing undesirable characters (eg, ":<>[]") with '_', avoid duplicates. Eg, "PaSR<psiReactionThermo>:Qdot" becomes "PaSR_psiReactionThermo_Qdot" instead of "PaSR_psiReactionThermo__Qdot" ENH: additional ensightCase::padded static method	2022-07-13 19:21:11 +02:00
				`@ -1,2 +0,0 @@`
				`auto parcels = parcelCloudModelList(g, rhoInf, U, muc);`
				`@ -1 +0,0 @@`
				`regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();`