WIP: electricPotential: add space charge density hook

ENH: ChargeCloud: new space charge density function object
ENH: ChargeCloud: add Coulomb force model
2023-02-24 08:34:08 +00:00 · 2023-02-24 08:34:08 +00:00 · 2023-02-24 08:34:08 +00:00 · 2023-02-24 08:34:08 +00:00 · 2023-02-24 08:34:06 +00:00 · 2023-02-24 08:33:08 +00:00
277 changed files with 7559 additions and 1761 deletions
--- a/applications/solvers/combustion/fireFoam/createClouds.H
+++ b/applications/solvers/combustion/fireFoam/createClouds.H
@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingCloud parcels
+reactingCloud parcels
 (
    "reactingCloud1",
+    g,
    rho,
    U,
-    g,
    slgThermo
 );
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@ -37,7 +37,7 @@ Description

 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicReactingCloud.H"
+#include "reactingCloud.H"
 #include "surfaceFilmModel.H"
 #include "pyrolysisModelCollection.H"
 #include "radiationModel.H"
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C
@ -42,11 +42,11 @@ Description
 #include "CorrectPhi.H"

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

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

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

        mesh.update();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        mesh.update();

        U.correctBoundaryConditions();

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

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

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

        runTime.write();

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

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

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

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

        laminarTransport.correct();

        mu = laminarTransport.nu()*rhoInfValue;

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

        runTime.write();

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

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

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

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

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

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

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

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

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

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

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

 #include "reactingParcelFoam.C"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        mesh.update();

        U.correctBoundaryConditions();

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

        runTime.write();

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

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

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

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

        runTime.write();

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

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

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

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

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

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

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

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

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

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

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

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

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

 surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
--- a/doc/Build.md
+++ b/doc/Build.md
@ -1,23 +1,23 @@
 ## Getting the code

-Links to all code packs are available on https://dl.openfoam.com. For OpenFOAM-v2206:
+Links to all code packs are available on https://dl.openfoam.com. For OpenFOAM-v2212:

 - https://dl.openfoam.com/source/latest/
- Source: https://dl.openfoam.com/source/v2206/OpenFOAM-v2206.tgz
- ThirdParty: https://dl.openfoam.com/source/v2206/ThirdParty-v2206.tgz
+- Source: https://dl.openfoam.com/source/v2212/OpenFOAM-v2212.tgz
+- ThirdParty: https://dl.openfoam.com/source/v2212/ThirdParty-v2212.tgz

 ## OpenFOAM&reg; Quick Build Guide

 Prior to building, ensure that the [system requirements][link openfoam-require]
 are satisfied (including any special [cross-compiling][wiki-cross-compile]
 considerations), and source the correct OpenFOAM environment.
-For example, for the OpenFOAM-v2206 version:
+For example, for the OpenFOAM-v2212 version:
 ```
-source <installation path>/OpenFOAM-v2206/etc/bashrc
+source <installation path>/OpenFOAM-v2212/etc/bashrc
 ```
 e.g. if installed under the `~/openfoam` directory
 ```
-source ~/openfoam/OpenFOAM-v2206/etc/bashrc
+source ~/openfoam/OpenFOAM-v2212/etc/bashrc
 ```


--- a/etc/bashrc
+++ b/etc/bashrc
@ -55,7 +55,7 @@

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

 #------------------------------------------------------------------------------
 # Configuration environment variables.
--- a/etc/cshrc
+++ b/etc/cshrc
@ -55,7 +55,7 @@

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

 #------------------------------------------------------------------------------
 # Configuration environment variables.
--- a/src/functionObjects/lagrangian/Make/options
+++ b/src/functionObjects/lagrangian/Make/options
@ -12,7 +12,14 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/faOptions/lnInclude \
+    \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude

 LIB_LIBS = \
    -lfiniteVolume \
--- a/src/functionObjects/lagrangian/cloudInfo/cloudInfo.C
+++ b/src/functionObjects/lagrangian/cloudInfo/cloudInfo.C
@ -28,7 +28,7 @@ License

 #include "cloudInfo.H"
 #include "cloud.H"
-#include "kinematicCloud.H"
+#include "parcelCloudModel.H"
 #include "dictionary.H"
 #include "mathematicalConstants.H"
 #include "PstreamReduceOps.H"
@ -147,7 +147,7 @@ bool Foam::functionObjects::cloudInfo::performAction(unsigned request)

        const word& cloudName = names()[cloudi];

-        const auto* kinCloudPtr = obr_.cfindObject<kinematicCloud>(cloudName);
+        const auto* kinCloudPtr = obr_.cfindObject<parcelCloudModel>(cloudName);

        if (!kinCloudPtr)
        {
--- a/src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.C
+++ b/src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.C
@ -92,10 +92,10 @@ Foam::functionObjects::icoUncoupledKinematicCloud::icoUncoupledKinematicCloud
    kinematicCloud_
    (
        kinematicCloudName_,
+        g_,
        rho_,
        U_,
-        mu_,
-        g_
+        mu_
    )
 {}

--- a/src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.H
+++ b/src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.H
@ -76,7 +76,7 @@ SourceFiles
 #define functionObjects_icoUncoupledKinematicCloud_H

 #include "fvMeshFunctionObject.H"
-#include "basicKinematicCollidingCloud.H"
+#include "kinematicCollidingCloud.H"
 #include "volFields.H"
 #include "uniformDimensionedFields.H"

@ -122,7 +122,7 @@ class icoUncoupledKinematicCloud
        word kinematicCloudName_;

        //- The kinematic cloud
-        basicKinematicCollidingCloud kinematicCloud_;
+        kinematicCollidingCloud kinematicCloud_;


    // Private member functions
--- a/src/functionObjects/solvers/electricPotential/electricPotential.C
+++ b/src/functionObjects/solvers/electricPotential/electricPotential.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -45,16 +45,20 @@ namespace functionObjects

 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

-Foam::volScalarField&
-Foam::functionObjects::electricPotential::operandField()
+Foam::volScalarField& Foam::functionObjects::electricPotential::getOrReadField
+(
+    const word& fieldName
+) const
 {
-    if (!foundObject<volScalarField>(fieldName_))
+    auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+    if (!ptr)
    {
-        auto tfldPtr = tmp<volScalarField>::New
+        ptr = new volScalarField
        (
            IOobject
            (
-                fieldName_,
+                fieldName,
                mesh_.time().timeName(),
                mesh_,
                IOobject::MUST_READ,
@ -62,10 +66,10 @@ Foam::functionObjects::electricPotential::operandField()
            ),
            mesh_
        );
-        store(fieldName_, tfldPtr);
+        mesh_.objectRegistry::store(ptr);
    }

-    return lookupObjectRef<volScalarField>(fieldName_);
+    return *ptr;
 }


@ -79,7 +83,7 @@ Foam::functionObjects::electricPotential::sigma() const
        mesh_.time(),
        IOobject::NO_READ,
        IOobject::NO_WRITE,
-        false
+        IOobject::NO_REGISTER
    );

    if (phases_.size())
@ -126,7 +130,7 @@ Foam::functionObjects::electricPotential::epsilonm() const
        mesh_.time(),
        IOobject::NO_READ,
        IOobject::NO_WRITE,
-        false
+        IOobject::NO_REGISTER
    );

    if (phases_.size())
@ -166,6 +170,7 @@ Foam::functionObjects::electricPotential::electricPotential
 )
 :
    fvMeshFunctionObject(name, runTime, dict),
+    fvOptions_(mesh_),
    phasesDict_(dict.subOrEmptyDict("phases")),
    phaseNames_(),
    phases_(),
@ -205,14 +210,16 @@ Foam::functionObjects::electricPotential::electricPotential
            IOobject::scopedName(typeName, "V")
        )
    ),
+    erhoName_(),
    nCorr_(1),
-    writeDerivedFields_(false)
+    writeDerivedFields_(false),
+    spaceChargeDensity_(false)
 {
    read(dict);

    // Force creation of transported field so any BCs using it can
    // look it up
-    volScalarField& eV = operandField();
+    volScalarField& eV = getOrReadField(fieldName_);
    eV.correctBoundaryConditions();
 }

@ -221,88 +228,101 @@ Foam::functionObjects::electricPotential::electricPotential

 bool Foam::functionObjects::electricPotential::read(const dictionary& dict)
 {
-    if (fvMeshFunctionObject::read(dict))
+    if (!fvMeshFunctionObject::read(dict))
    {
-        Log << type() << " read: " << name() << endl;
+        return false;
+    }

-        dict.readIfPresent("sigma", sigma_);
-        dict.readIfPresent("epsilonr", epsilonr_);
-        dict.readIfPresent("nCorr", nCorr_);
-        dict.readIfPresent("writeDerivedFields", writeDerivedFields_);
+    Log << type() << " read: " << name() << endl;

-        // If flow is multiphase
-        if (!phasesDict_.empty())
+    {
+        const dictionary* dictptr = dict.findDict("fvOptions");
+        if (dictptr)
        {
-            phaseNames_.setSize(phasesDict_.size());
-            phases_.setSize(phasesDict_.size());
-            sigmas_.setSize(phasesDict_.size());
+            fvOptions_.reset(*dictptr);
+        }
+    }
+
+    dict.readIfPresent("sigma", sigma_);
+    dict.readIfPresent("epsilonr", epsilonr_);
+    dict.readIfPresent("nCorr", nCorr_);
+    dict.readIfPresent("writeDerivedFields", writeDerivedFields_);
+
+    if (dict.readIfPresent("spaceChargeDensity", spaceChargeDensity_))
+    {
+        dict.readEntry("erho", erhoName_);
+    }
+
+    // If flow is multiphase
+    if (!phasesDict_.empty())
+    {
+        phaseNames_.setSize(phasesDict_.size());
+        phases_.setSize(phasesDict_.size());
+        sigmas_.setSize(phasesDict_.size());
+
+        if (writeDerivedFields_)
+        {
+            epsilonrs_.setSize(phasesDict_.size());
+        }
+
+        label phasei = 0;
+        for (const entry& dEntry : phasesDict_)
+        {
+            const word& key = dEntry.keyword();
+
+            if (!dEntry.isDict())
+            {
+                FatalIOErrorInFunction(phasesDict_)
+                    << "Entry " << key << " is not a dictionary" << nl
+                    << exit(FatalIOError);
+            }
+
+            const dictionary& subDict = dEntry.dict();
+
+            phaseNames_[phasei] = key;
+
+            sigmas_.set
+            (
+                phasei,
+                new dimensionedScalar
+                (
+                    sqr(dimCurrent)*pow3(dimTime)/(dimMass*pow3(dimLength)),
+                    subDict.getCheck<scalar>
+                    (
+                        "sigma",
+                        scalarMinMax::ge(SMALL)
+                    )
+                )
+            );

            if (writeDerivedFields_)
            {
-                epsilonrs_.setSize(phasesDict_.size());
-            }
-
-            label phasei = 0;
-            for (const entry& dEntry : phasesDict_)
-            {
-                const word& key = dEntry.keyword();
-
-                if (!dEntry.isDict())
-                {
-                    FatalIOErrorInFunction(phasesDict_)
-                        << "Entry " << key << " is not a dictionary" << nl
-                        << exit(FatalIOError);
-                }
-
-                const dictionary& subDict = dEntry.dict();
-
-                phaseNames_[phasei] = key;
-
-                sigmas_.set
+                epsilonrs_.set
                (
                    phasei,
                    new dimensionedScalar
                    (
-                        sqr(dimCurrent)*pow3(dimTime)/(dimMass*pow3(dimLength)),
+                        dimless,
                        subDict.getCheck<scalar>
                        (
-                            "sigma",
+                            "epsilonr",
                            scalarMinMax::ge(SMALL)
                        )
                    )
                );
-
-                if (writeDerivedFields_)
-                {
-                    epsilonrs_.set
-                    (
-                        phasei,
-                        new dimensionedScalar
-                        (
-                            dimless,
-                            subDict.getCheck<scalar>
-                            (
-                                "epsilonr",
-                                scalarMinMax::ge(SMALL)
-                            )
-                        )
-                    );
-                }
-
-                ++phasei;
            }

-            forAll(phaseNames_, i)
-            {
-                phases_.set
-                (
-                    i,
-                    mesh_.getObjectPtr<volScalarField>(phaseNames_[i])
-                );
-            }
+            ++phasei;
        }

-        return true;
+        forAll(phaseNames_, i)
+        {
+            phases_.set
+            (
+                i,
+                mesh_.getObjectPtr<volScalarField>(phaseNames_[i])
+            );
+        }
    }

    return false;
@ -316,17 +336,31 @@ bool Foam::functionObjects::electricPotential::execute()
    tmp<volScalarField> tsigma = this->sigma();
    const volScalarField& sigma = tsigma();

-    volScalarField& eV = operandField();
+    tmp<volScalarField> tepsilon = this->epsilonm();
+    const auto& epsilon = tepsilon();
+
+    volScalarField& eV = getOrReadField(fieldName_);

    for (label i = 1; i <= nCorr_; ++i)
    {
        fvScalarMatrix eVEqn
        (
          - fvm::laplacian(sigma, eV)
+         ==
+            fvOptions_(eV)
        );

+        if (spaceChargeDensity_)
+        {
+            volScalarField& erho = getOrReadField(erhoName_);
+
+            eVEqn += erho/epsilon;
+        }
+
        eVEqn.relax();

+        fvOptions_.constrain(eVEqn);
+
        eVEqn.solve();
    }

@ -342,7 +376,7 @@ bool Foam::functionObjects::electricPotential::write()
        << tab << fieldName_
        << endl;

-    volScalarField& eV = operandField();
+    volScalarField& eV = getOrReadField(fieldName_);

    if (writeDerivedFields_)
    {
@ -356,7 +390,7 @@ bool Foam::functionObjects::electricPotential::write()
                mesh_.time(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
-                false
+                IOobject::NO_REGISTER
            ),
            -fvc::grad(eV),
            calculatedFvPatchField<vector>::typeName
@ -379,7 +413,7 @@ bool Foam::functionObjects::electricPotential::write()
                mesh_.time(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
-                false
+                IOobject::NO_REGISTER
            ),
            -tsigma*fvc::grad(eV),
            calculatedFvPatchField<vector>::typeName
@ -402,7 +436,7 @@ bool Foam::functionObjects::electricPotential::write()
                mesh_.time(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
-                false
+                IOobject::NO_REGISTER
            ),
            fvc::div(tepsilonm*E),
            calculatedFvPatchField<scalar>::typeName
--- a/src/functionObjects/solvers/electricPotential/electricPotential.H
+++ b/src/functionObjects/solvers/electricPotential/electricPotential.H
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -37,7 +37,7 @@ Description
    The steady-state equation of the charge conservation:

    \f[
-        \nabla \cdot \left( \sigma \nabla V \right) = 0
+        \nabla \cdot \left(\sigma \nabla V \right) = 0
    \f]

    where
@ -121,6 +121,8 @@ Usage
        nCorr                 <label>;
        writeDerivedFields    <bool>;
        fieldName             <word>;
+        spaceChargeDensity    <bool>;
+        erhoName              <word>;

        // Inherited entries
        ...
@ -137,6 +139,9 @@ Usage
      nCorr     | Number of corrector iterations     | label | no   | 1
      writeDerivedFields | Flag to write extra fields | bool | no   | false
      fieldName | Name of operand field   | word | no | electricPotential:V
+      spaceChargeDensity | Flag to include space charge density | bool | no <!--
+                -->                                                 | false
+      erhoName  | Name of space charge density field  | word | no   | -
    \endtable

    The inherited entries are elaborated in:
@ -160,6 +165,7 @@ SourceFiles

 #include "fvMeshFunctionObject.H"
 #include "volFields.H"
+#include "fvOptionList.H"

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

@ -178,6 +184,9 @@ class electricPotential
 {
    // Private Data

+        //- Run-time selectable finite volume options, e.g. sources, constraints
+        fv::optionList fvOptions_;
+
        //- Dictionary of phase data
        dictionary phasesDict_;

@ -202,6 +211,9 @@ class electricPotential
        //- Name of the operand field
        word fieldName_;

+        //- Name of the space charge density field
+        word erhoName_;
+
        //- Number of corrector iterations
        label nCorr_;

@ -209,11 +221,16 @@ class electricPotential
        //- electric field, current density and free-charge density
        bool writeDerivedFields_;

+        //- Flag to add the space charge density to the rhs of the potential eq
+        bool spaceChargeDensity_;
+

    // Private Member Functions

-        //- Return reference to the registered operand field
-        volScalarField& operandField();
+        //- Return requested field from the object registry
+        //- or read+register the field to the object registry
+        volScalarField& getOrReadField(const word& fieldName) const;
+

        //- Return the isotropic electrical conductivity field of the mixture
        tmp<volScalarField> sigma() const;
--- a/src/lagrangian/Allwmake
+++ b/src/lagrangian/Allwmake
@ -9,9 +9,7 @@ wmake $targetType basic
 wmake $targetType solidParticle
 wmake $targetType intermediate
 wmake $targetType turbulence
-wmake $targetType spray
 wmake $targetType DSMC
-wmake $targetType coalCombustion

 molecularDynamics/Allwmake $targetType $*

--- a/src/lagrangian/coalCombustion/Make/files
+++ b/src/lagrangian/coalCombustion/Make/files
@ -1,5 +0,0 @@
-/* Coal parcel and sub-models */
-coalParcel/makeCoalParcelSubmodels.C
-coalCloudList/coalCloudList.C
-
-LIB = $(FOAM_LIBBIN)/libcoalCombustion
--- a/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
+++ b/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
@ -1,264 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
-------------------------------------------------------------------------------
-    Copyright (C) 2012-2016 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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvMatrices.H"
-#include "volFields.H"
-#include "DimensionedField.H"
-
-Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh>>
-Foam::coalCloudList::UTrans() const
-{
-    tmp<volVectorField::Internal> tfld
-    (
-        new volVectorField::Internal
-        (
-            IOobject
-            (
-                "UTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedVector(dimMass*dimVelocity, Zero)
-        )
-    );
-
-    volVectorField::Internal& fld = tfld.ref();
-
-    forAll(*this, i)
-    {
-        fld += operator[](i).UTrans();
-    }
-
-    return tfld;
-}
-
-
-Foam::tmp<Foam::fvVectorMatrix> Foam::coalCloudList::SU
-(
-    volVectorField& U
-) const
-{
-    tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
-    fvVectorMatrix& fvm = tfvm.ref();
-
-    forAll(*this, i)
-    {
-        fvm += operator[](i).SU(U);
-    }
-
-    return tfvm;
-}
-
-
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::coalCloudList::hsTrans() const
-{
-    tmp<volScalarField::Internal> tfld
-    (
-        new volScalarField::Internal
-        (
-            IOobject
-            (
-                "hsTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(dimEnergy, Zero)
-        )
-    );
-
-    volScalarField::Internal& fld = tfld.ref();
-
-    forAll(*this, i)
-    {
-        fld += operator[](i).hsTrans();
-    }
-
-    return tfld;
-}
-
-
-Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Sh
-(
-    volScalarField& hs
-) const
-{
-    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
-    fvScalarMatrix& fvm = tfvm.ref();
-
-    forAll(*this, i)
-    {
-        fvm += operator[](i).Sh(hs);
-    }
-
-    return tfvm;
-}
-
-
-Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::SYi
-(
-    const label ii,
-    volScalarField& Yi
-) const
-{
-    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
-    fvScalarMatrix& fvm = tfvm.ref();
-
-    forAll(*this, i)
-    {
-        fvm += operator[](i).SYi(ii, Yi);
-    }
-
-    return tfvm;
-}
-
-
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::coalCloudList::rhoTrans() const
-{
-    tmp<volScalarField::Internal> tfld
-    (
-        new volScalarField::Internal
-        (
-            IOobject
-            (
-                "rhoTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(dimMass, Zero)
-        )
-    );
-
-    volScalarField::Internal& fld = tfld.ref();
-
-    forAll(*this, i)
-    {
-        forAll(operator[](i).rhoTrans(), j)
-        {
-            fld += operator[](i).rhoTrans()[j];
-        }
-    }
-
-    return tfld;
-}
-
-
-
-
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::coalCloudList::Srho() const
-{
-    tmp<volScalarField::Internal> tfld
-    (
-        new volScalarField::Internal
-        (
-            IOobject
-            (
-                "rhoTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(dimDensity/dimTime, Zero)
-        )
-    );
-
-    volScalarField::Internal& fld = tfld.ref();
-
-    forAll(*this, i)
-    {
-        fld += operator[](i).Srho();
-    }
-
-    return tfld;
-}
-
-
-Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
-Foam::coalCloudList::Srho
-(
-    const label i
-) const
-{
-    tmp<volScalarField::Internal> tfld
-    (
-        new volScalarField::Internal
-        (
-            IOobject
-            (
-                "rhoTransEff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(dimDensity/dimTime, Zero)
-        )
-    );
-
-    volScalarField::Internal& fld = tfld.ref();
-
-    forAll(*this, j)
-    {
-        fld += operator[](j).Srho(i);
-    }
-
-    return tfld;
-}
-
-
-Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Srho
-(
-    volScalarField& rho
-) const
-{
-    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(rho, dimMass/dimTime));
-    fvScalarMatrix& fvm = tfvm.ref();
-
-    forAll(*this, i)
-    {
-        fvm += operator[](i).Srho(rho);
-    }
-
-    return tfvm;
-}
-
-
-// ************************************************************************* //
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@ -2,59 +2,81 @@ PARCELS=parcels
 BASEPARCELS=$(PARCELS)/baseClasses
 DERIVEDPARCELS=$(PARCELS)/derived

+
 CLOUDS=clouds
-BASECLOUDS=$(CLOUDS)/baseClasses
+TEMPLATECLOUDS=$(CLOUDS)/Templates
+$(TEMPLATECLOUDS)/KinematicCloud/KinematicCloudName.C
+$(TEMPLATECLOUDS)/ThermoCloud/ThermoCloudName.C
+$(TEMPLATECLOUDS)/ReactingCloud/ReactingCloudName.C
+$(TEMPLATECLOUDS)/ReactingMultiphaseCloud/ReactingMultiphaseCloudName.C
+$(TEMPLATECLOUDS)/ChargeCloud/ChargeCloudName.C
+$(TEMPLATECLOUDS)/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudName.C
+$(TEMPLATECLOUDS)/CollidingCloud/CollidingCloudName.C
+$(TEMPLATECLOUDS)/MPPICCloud/MPPICCloudName.C
+$(TEMPLATECLOUDS)/SprayCloud/SprayCloudName.C
+
+
 DERIVEDCLOUDS=$(CLOUDS)/derived
-
-
-/* Cloud base classes */
-$(BASECLOUDS)/kinematicCloud/kinematicCloud.C
-$(BASECLOUDS)/thermoCloud/thermoCloud.C
-$(BASECLOUDS)/reactingCloud/reactingCloud.C
-$(BASECLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
-$(BASECLOUDS)/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C
+$(DERIVEDCLOUDS)/kinematicCloud/kinematicCloud.C
+$(DERIVEDCLOUDS)/thermoCloud/thermoCloud.C
+$(DERIVEDCLOUDS)/reactingCloud/reactingCloud.C
+$(DERIVEDCLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
+$(DERIVEDCLOUDS)/chargeCloud/chargeCloud.C
+$(DERIVEDCLOUDS)/heterogeneousReactingCloud/heterogeneousReactingCloud.C
+$(DERIVEDCLOUDS)/kinematicCollidingCloud/kinematicCollidingCloud.C
+$(DERIVEDCLOUDS)/kinematicMPPICCloud/kinematicMPPICCloud.C
+$(DERIVEDCLOUDS)/sprayCloud/sprayCloud.C


 /* kinematic parcel sub-models */
-KINEMATICPARCEL=$(DERIVEDPARCELS)/basicKinematicParcel
-$(KINEMATICPARCEL)/defineBasicKinematicParcel.C
-$(KINEMATICPARCEL)/makeBasicKinematicParcelSubmodels.C
+KINEMATICPARCEL=$(DERIVEDPARCELS)/kinematicParcel
+$(KINEMATICPARCEL)/defineKinematicParcel.C
+$(KINEMATICPARCEL)/makeKinematicParcelSubmodels.C


 /* kinematic colliding parcel sub-models */
-KINEMATICCOLLIDINGPARCEL=$(DERIVEDPARCELS)/basicKinematicCollidingParcel
-$(KINEMATICCOLLIDINGPARCEL)/defineBasicKinematicCollidingParcel.C
-$(KINEMATICCOLLIDINGPARCEL)/makeBasicKinematicCollidingParcelSubmodels.C
+KINEMATICCOLLIDINGPARCEL=$(DERIVEDPARCELS)/kinematicCollidingParcel
+$(KINEMATICCOLLIDINGPARCEL)/defineKinematicCollidingParcel.C
+$(KINEMATICCOLLIDINGPARCEL)/makeKinematicCollidingParcelSubmodels.C


 /* thermo parcel sub-models */
-THERMOPARCEL=$(DERIVEDPARCELS)/basicThermoParcel
-$(THERMOPARCEL)/defineBasicThermoParcel.C
-$(THERMOPARCEL)/makeBasicThermoParcelSubmodels.C
+THERMOPARCEL=$(DERIVEDPARCELS)/thermoParcel
+$(THERMOPARCEL)/defineThermoParcel.C
+$(THERMOPARCEL)/makeThermoParcelSubmodels.C


 /* reacting parcel sub-models */
-REACTINGPARCEL=$(DERIVEDPARCELS)/basicReactingParcel
-$(REACTINGPARCEL)/defineBasicReactingParcel.C
-$(REACTINGPARCEL)/makeBasicReactingParcelSubmodels.C
+REACTINGPARCEL=$(DERIVEDPARCELS)/reactingParcel
+$(REACTINGPARCEL)/defineReactingParcel.C
+$(REACTINGPARCEL)/makeReactingParcelSubmodels.C


 /* reacting multiphase parcel sub-models */
-REACTINGMPPARCEL=$(DERIVEDPARCELS)/basicReactingMultiphaseParcel
-$(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
-$(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
+REACTINGMPPARCEL=$(DERIVEDPARCELS)/reactingMultiphaseParcel
+$(REACTINGMPPARCEL)/defineReactingMultiphaseParcel.C
+$(REACTINGMPPARCEL)/makeReactingMultiphaseParcelSubmodels.C


-/* heterogeous reacting parcel sub-models */
-REACTINGHETERMPPARCEL=$(DERIVEDPARCELS)/basicHeterogeneousReactingParcel
-$(REACTINGHETERMPPARCEL)/defineBasicHeterogeneousReactingParcel.C
-$(REACTINGHETERMPPARCEL)/makeBasicHeterogeneousReactingParcelSubmodels.C
+/* charge parcel sub-models */
+CHARGEPARCEL=$(DERIVEDPARCELS)/chargeParcel
+$(CHARGEPARCEL)/defineChargeParcel.C
+$(CHARGEPARCEL)/makeChargeParcelSubmodels.C

+/* heterogeneous reacting parcel sub-models */
+REACTINGHETERMPPARCEL=$(DERIVEDPARCELS)/heterogeneousReactingParcel
+$(REACTINGHETERMPPARCEL)/defineHeterogeneousReactingParcel.C
+$(REACTINGHETERMPPARCEL)/makeHeterogeneousReactingParcelSubmodels.C

 /* kinematic MPPIC parcel sub-models */
-KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/basicKinematicMPPICParcel
-$(KINEMATICMPPICPARCEL)/defineBasicKinematicMPPICParcel.C
-$(KINEMATICMPPICPARCEL)/makeBasicKinematicMPPICParcelSubmodels.C
+KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/kinematicMPPICParcel
+$(KINEMATICMPPICPARCEL)/defineKinematicMPPICParcel.C
+$(KINEMATICMPPICPARCEL)/makeKinematicMPPICParcelSubmodels.C
+
+/* spray parcel sub-models */
+SPRAYPARCEL=$(DERIVEDPARCELS)/sprayParcel
+$(SPRAYPARCEL)/defineSprayParcel.C
+$(SPRAYPARCEL)/makeSprayParcelSubmodels.C


 /* bolt-on models */
@ -86,6 +108,11 @@ $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphase
 $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphaseParcelInjectionDataIO.C
 $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphaseParcelInjectionDataIOList.C

+CHARGEINJECTION=submodels/Charge/InjectionModel
+$(CHARGEINJECTION)/ChargeLookupTableInjection/chargeParcelInjectionData.C
+$(CHARGEINJECTION)/ChargeLookupTableInjection/chargeParcelInjectionDataIO.C
+$(CHARGEINJECTION)/ChargeLookupTableInjection/chargeParcelInjectionDataIOList.C
+
 MPPICPARTICLESTRESS=submodels/MPPIC/ParticleStressModels
 $(MPPICPARTICLESTRESS)/ParticleStressModel/ParticleStressModel.C
 $(MPPICPARTICLESTRESS)/HarrisCrighton/HarrisCrighton.C
@ -130,4 +157,10 @@ conversion/ensight/ensightOutputCloud.C
 conversion/vtk/foamVtkLagrangianWriter.C


+clouds/baseCloudInterface.C
+clouds/parcelCloudModel.C
+clouds/parcelCloudModelNew.C
+clouds/parcelCloudModelList/parcelCloudModelList.C
+
+
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate
--- a/src/lagrangian/intermediate/clouds/BaseCloud.H
+++ b/src/lagrangian/intermediate/clouds/BaseCloud.H
@ -0,0 +1,275 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2021 OpenFOAM Foundation
+    Copyright (C) 2022 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/>.
+
+Class
+    Foam::BaseCloud
+
+Description
+    Base template for parcel clouds. Inserts the baseCloudInterface
+    virtualisation layer into the class. Also defines default zero-return
+    source methods to enable the less functional clouds to be used in more
+    complex situations.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_BaseCloud_H
+#define Foam_BaseCloud_H
+
+#include "Cloud.H"
+#include "baseCloudInterface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class SLGThermo;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+class BaseCloud
+:
+    public Cloud<ParticleType>,
+    virtual public baseCloudInterface
+{
+    // Private data
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+
+public:
+
+    // Constructors
+
+        BaseCloud
+        (
+            const word& cloudName,
+            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const volScalarField& mu,
+            const bool readFields = true
+        )
+        :
+            Cloud<ParticleType>(rho.mesh(), cloudName, false),
+            mesh_(rho.mesh())
+        {}
+
+
+        BaseCloud
+        (
+            const word& cloudName,
+            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const SLGThermo& thermo,
+            const bool readFields = true
+        )
+        :
+            Cloud<ParticleType>(rho.mesh(), cloudName, false),
+            mesh_(rho.mesh())
+        {}
+
+        //- Copy constructor with new name
+        BaseCloud
+        (
+            BaseCloud<ParticleType>& c,
+            const word& name
+        )
+        :
+            Cloud<ParticleType>(c.mesh(), name, c),
+            mesh_(c.mesh_)
+        {}
+
+        //- Copy constructor with new name - creates bare cloud
+        BaseCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const BaseCloud<ParticleType>& c
+        )
+        :
+            Cloud<ParticleType>(mesh, name, IDLList<ParticleType>()),
+            mesh_(mesh)
+        {}
+
+    //- Destructor
+    virtual ~BaseCloud() = default;
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return reference to the mesh
+            const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+
+
+        // Sources
+
+            // Momentum
+
+                //- Momentum transfer [kg m/s]
+                virtual tmp<volVectorField::Internal> UTrans() const
+                {
+                    return volVectorField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), ":UTrans"),
+                        mesh_,
+                        dimensionedVector(dimMass*dimVelocity, Zero)
+                    );
+                }
+
+                //- Momentum transfer coefficient [kg]
+                virtual tmp<volScalarField::Internal> UCoeff() const
+                {
+                    return volScalarField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), "UCoeff"),
+                        mesh_,
+                        dimensionedScalar(dimMass, Zero)
+                    );
+                }
+
+                //- Return momentum source term [kg m/s^2]
+                virtual tmp<fvVectorMatrix> SU(volVectorField& U) const
+                {
+                    return tmp<fvVectorMatrix>::New(U, dimMass*dimAcceleration);
+                }
+
+
+            // Energy
+
+                //- Return sensible enthalpy transfer [J]
+                virtual tmp<volScalarField::Internal> hsTrans() const
+                {
+                    return volScalarField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), "hsTrans"),
+                        mesh_,
+                        dimensionedScalar(dimEnergy, Zero)
+                    );
+                }
+
+                //- Return sensible enthalpy transfer coefficient [J/K]
+                virtual tmp<volScalarField::Internal> hsCoeff() const
+                {
+                    return volScalarField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), "hsCoeff"),
+                        mesh_,
+                        dimensionedScalar(dimEnergy/dimTemperature, Zero)
+                    );
+                }
+
+                //- Return sensible enthalpy source term [J/s]
+                virtual tmp<fvScalarMatrix> Sh(volScalarField& hs) const
+                {
+                    return tmp<fvScalarMatrix>::New(hs, dimEnergy/dimTime);
+                }
+
+
+            // Radiation - overrides thermoCloud virtual abstract members
+
+                //- Return equivalent particulate emission [kg/m/s^3]
+                virtual tmp<volScalarField> Ep() const
+                {
+                    return volScalarField::New
+                    (
+                        IOobject::scopedName(this->name(), "radiation_Ep"),
+                        mesh_,
+                        dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
+                    );
+                }
+
+                //- Return equivalent particulate absorption [1/m]
+                virtual tmp<volScalarField> ap() const
+                {
+                    return volScalarField::New
+                    (
+                        IOobject::scopedName(this->name(), "radiation_ap"),
+                        mesh_,
+                        dimensionedScalar(dimless/dimLength, Zero)
+                    );
+                }
+
+                //- Return equivalent particulate scattering factor [1/m]
+                virtual tmp<volScalarField> sigmap() const
+                {
+                    return volScalarField::New
+                    (
+                        IOobject::scopedName(this->name(), "radiation_sigmap"),
+                        mesh_,
+                        dimensionedScalar(dimless/dimLength, Zero)
+                    );
+                }
+
+
+            // Mass
+
+                //- Return mass source term for specie [kg/s]
+                virtual tmp<fvScalarMatrix> SYi
+                (
+                    const label speciei,
+                    volScalarField& Yi
+                ) const
+                {
+                    return tmp<fvScalarMatrix>::New(Yi, dimMass/dimTime);
+                }
+
+                //- Return total mass source term [kg/m^3/s]
+                virtual tmp<volScalarField::Internal> Srho() const
+                {
+                    return volScalarField::Internal::New
+                    (
+                        IOobject::scopedName(this->name(), "Srho"),
+                        mesh_,
+                        dimensionedScalar(dimDensity/dimTime, Zero)
+                    );
+                }
+
+                //- Return total mass source term [kg/s]
+                virtual tmp<fvScalarMatrix> Srho(volScalarField& rho) const
+                {
+                    return tmp<fvScalarMatrix>::New(rho, dimMass/dimTime);
+                }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/baseClasses/thermoCloud/thermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/baseClasses/thermoCloud/thermoCloud.H
@ -5,7 +5,8 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2020-2021 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -24,65 +25,67 @@ License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

 Class
-    Foam::thermoCloud
+    Foam::ParcelCloudModel

 Description
-    Virtual abstract base class for templated ThermoCloud
+    Outermost template for parcel clouds. Adds the parcelCloudModel
+    virtualisation layer and forwards the methods required by that layer.

 SourceFiles
-    thermoCloud.C
+    ParcelCloudModel.C

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

-#ifndef thermoCloud_H
-#define thermoCloud_H
+#ifndef Foam_ParcelCloudModel_H
+#define Foam_ParcelCloudModel_H

-#include "volFields.H"
+#include "parcelCloudModel.H"

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

 namespace Foam
 {

-/*---------------------------------------------------------------------------*\
-                       Class thermoCloud Declaration
-\*---------------------------------------------------------------------------*/
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-class thermoCloud
+template<class CloudType>
+class ParcelCloudModel
+:
+    public parcelCloudModel,
+    public CloudType
 {
 public:

-    //- Runtime type information
-    TypeName("thermoCloud");
-
-
    // Constructors

-        //- Null constructor
-        thermoCloud() = default;
-
-
-    //- Destructor
-    virtual ~thermoCloud() = default;
+        //- Inherit base constructors
+        using CloudType::CloudType;


    // Member Functions

-        // Access
+        //- Evolve the cloud
+        virtual void evolveME() final
+        {
+            // CloudType::evolve();
+            if (this->solution().canEvolve())
+            {
+                typedef typename CloudType::particleType parcelType;

-            // Radiation
+                typename parcelType::trackingData td(*this);

-                //- Return tmp equivalent particulate emission
-                virtual tmp<volScalarField> Ep() const = 0;
+                this->solve(*this, td);
+            }
+        }

-                //- Return tmp equivalent particulate absorption
-                virtual tmp<volScalarField> ap() const = 0;
-
-                //- Return tmp equivalent particulate scattering factor
-                virtual tmp<volScalarField> sigmap() const = 0;
+        //- Call this before a topology change. Stores the particles global
+        //- positions in the database for use during mapping.
+        virtual void storeGlobalPositionsME() const final
+        {
+            CloudType::storeGlobalPositions();
+        }
 };

-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 } // End namespace Foam
--- a/src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloud.C
@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ChargeCloud.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::setModels()
+{
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::cloudReset
+(
+    ChargeCloud<CloudType>& c
+)
+{
+    CloudType::cloudReset(c);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ChargeCloud<CloudType>::ChargeCloud
+(
+    const word& cloudName,
+    const dimensionedVector& g,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    CloudType(cloudName, g, rho, U, thermo, false),
+    cloudCopyPtr_(nullptr),
+    constProps_(this->particleProperties())
+{
+    if (this->solution().active())
+    {
+        setModels();
+
+        if (readFields)
+        {
+            parcelType::readFields(*this, this->composition());
+            this->deleteLostParticles();
+        }
+    }
+
+    if (this->solution().resetSourcesOnStartup())
+    {
+        resetSourceTerms();
+    }
+}
+
+
+template<class CloudType>
+Foam::ChargeCloud<CloudType>::ChargeCloud
+(
+    ChargeCloud<CloudType>& c,
+    const word& name
+)
+:
+    CloudType(c, name),
+    cloudCopyPtr_(nullptr),
+    constProps_(c.constProps_)
+{}
+
+
+template<class CloudType>
+Foam::ChargeCloud<CloudType>::ChargeCloud
+(
+    const fvMesh& mesh,
+    const word& name,
+    const ChargeCloud<CloudType>& c
+)
+:
+    CloudType(mesh, name, c),
+    cloudCopyPtr_(nullptr),
+    constProps_()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::checkParcelProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+
+    if (fullyDescribed)
+    {
+
+    }
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::storeState()
+{
+    cloudCopyPtr_.reset
+    (
+        static_cast<ChargeCloud<CloudType>*>
+        (
+            clone(this->name() + "Copy").ptr()
+        )
+    );
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::restoreState()
+{
+    cloudReset(cloudCopyPtr_());
+    cloudCopyPtr_.clear();
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::resetSourceTerms()
+{
+    CloudType::resetSourceTerms();
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::evolve()
+{
+    if (this->solution().canEvolve())
+    {
+        typename parcelType::trackingData td(*this);
+
+        this->solve(*this, td);
+    }
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::autoMap
+(
+    const mapPolyMesh& mapper
+)
+{
+    Cloud<parcelType>::autoMap(mapper);
+
+    this->updateMesh();
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::info()
+{
+    CloudType::info();
+}
+
+/*
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::writeFields() const
+{
+    if (this->compositionModel_)
+    {
+        CloudType::particleType::writeFields(*this, this->composition());
+    }
+}
+*/
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloud.H
@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+Class
+    Foam::ChargeCloud
+
+Group
+    grpLagrangianIntermediateClouds
+
+Description
+    Templated base class for electric-charge cloud
+
+    - Adds to multiphase reacting cloud
+      - electric charge
+
+SourceFiles
+    ChargeCloudI.H
+    ChargeCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_ChargeCloud_H
+#define Foam_ChargeCloud_H
+
+#include "className.H"
+#include "autoPtr.H"
+#include "Cloud.H"
+#include "volFields.H"
+#include "fvMatricesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+class SLGThermo;
+
+TemplateName(ChargeCloud);
+
+/*---------------------------------------------------------------------------*\
+                        Class ChargeCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ChargeCloud
+:
+    public CloudType,
+    public ChargeCloudName
+{
+public:
+
+    // Public Typedefs
+
+        //- Type of cloud this cloud was instantiated for
+        typedef CloudType cloudType;
+
+        //- Type of parcel the cloud was instantiated for
+        typedef typename CloudType::particleType parcelType;
+
+        //- Convenience typedef for this cloud type
+        typedef ChargeCloud<CloudType> chargeCloudType;
+
+
+private:
+
+    // Private Data
+
+        //- Cloud copy pointer
+        autoPtr<ChargeCloud<CloudType>> cloudCopyPtr_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        ChargeCloud(const ChargeCloud&) = delete;
+
+        //- No copy assignment
+        void operator=(const ChargeCloud&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        //- Parcel constant properties
+        typename parcelType::constantProperties constProps_;
+
+
+    // Protected Member Functions
+
+        // Initialisation
+
+            //- Set cloud sub-models
+            void setModels();
+
+
+        // Cloud Evolution Functions
+
+            //- Reset state of cloud
+            void cloudReset(ChargeCloud<CloudType>& c);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given carrier gas fields
+        ChargeCloud
+        (
+            const word& cloudName,
+            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const SLGThermo& thermo,
+            bool readFields = true
+        );
+
+        //- Copy constructor with new name
+        ChargeCloud
+        (
+            ChargeCloud<CloudType>& c,
+            const word& name
+        );
+
+        //- Copy constructor with new name - creates bare cloud
+        ChargeCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const ChargeCloud<CloudType>& c
+        );
+
+        //- Construct and return clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType>> clone(const word& name)
+        {
+            return autoPtr<Cloud<parcelType>>
+            (
+                new ChargeCloud(*this, name)
+            );
+        }
+
+        //- Construct and return bare clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType>> cloneBare(const word& name) const
+        {
+            return autoPtr<Cloud<parcelType>>
+            (
+                new ChargeCloud(this->mesh(), name, *this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ChargeCloud() = default;
+
+
+    // Member Functions
+
+        // Access
+
+        //- Return a reference to the cloud copy
+        inline const ChargeCloud& cloudCopy() const;
+
+        //- Return the constant properties
+        inline const typename parcelType::constantProperties&
+            constProps() const;
+
+        //- Return access to the constant properties
+        inline typename parcelType::constantProperties& constProps();
+
+
+        // Sub-models
+
+
+        // Cloud Evolution
+
+        //- Set parcel thermo properties
+        void setParcelThermoProperties
+        (
+            parcelType& parcel,
+            const scalar lagrangianDt
+        );
+
+        //- Check parcel properties
+        void checkParcelProperties
+        (
+            parcelType& parcel,
+            const scalar lagrangianDt,
+            const bool fullyDescribed
+        );
+
+        //- Store the current cloud state
+        void storeState();
+
+        //- Reset the current cloud to the previously stored state
+        void restoreState();
+
+        //- Reset the cloud source terms
+        void resetSourceTerms();
+
+        //- Evolve the cloud
+        void evolve();
+
+
+        // Mapping
+
+        //- Remap the cells of particles corresponding to the
+        //- mesh topology change with a default tracking data object
+        virtual void autoMap(const mapPolyMesh&);
+
+
+        // I-O
+
+        //- Print cloud information
+        void info();
+
+        //- Write the field data for the cloud
+        // virtual void writeFields() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ChargeCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ChargeCloud.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloudI.H
@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+inline const Foam::ChargeCloud<CloudType>&
+Foam::ChargeCloud<CloudType>::cloudCopy() const
+{
+    return *cloudCopyPtr_;
+}
+
+
+template<class CloudType>
+inline const typename CloudType::particleType::constantProperties&
+Foam::ChargeCloud<CloudType>::constProps() const
+{
+    return constProps_;
+}
+
+
+template<class CloudType>
+inline typename CloudType::particleType::constantProperties&
+Foam::ChargeCloud<CloudType>::constProps()
+{
+    return constProps_;
+}
+
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C
+++ b/src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -25,14 +25,13 @@ License

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

-#include "reactingHeterogeneousCloud.H"
+#include "ChargeCloud.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

 namespace Foam
 {
-    defineTypeNameAndDebug(reactingHeterogeneousCloud, 0);
+defineTypeNameAndDebug(ChargeCloudName, 0);
 }

-
 // ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C
@ -28,6 +28,7 @@ License
 #include "CollidingCloud.H"
 #include "CollisionModel.H"
 #include "NoCollision.H"
+#include "subCycleTime.H"

 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //

@ -88,14 +89,14 @@ template<class CloudType>
 Foam::CollidingCloud<CloudType>::CollidingCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
    const volScalarField& mu,
-    const dimensionedVector& g,
    bool readFields
 )
 :
-    CloudType(cloudName, rho, U, mu, g, false),
+    CloudType(cloudName, g, rho, U, mu, false),
    constProps_(this->particleProperties()),
    collisionModel_(nullptr)
 {
@ -148,13 +149,6 @@ Foam::CollidingCloud<CloudType>::CollidingCloud
 {}


-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::CollidingCloud<CloudType>::~CollidingCloud()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

 template<class CloudType>
--- a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H
@ -58,6 +58,8 @@ namespace Foam
 template<class CloudType>
 class CollisionModel;

+TemplateName(CollidingCloud);
+
 /*---------------------------------------------------------------------------*\
                       Class CollidingCloud Declaration
 \*---------------------------------------------------------------------------*/
@ -65,7 +67,8 @@ class CollisionModel;
 template<class CloudType>
 class CollidingCloud
 :
-    public CloudType
+    public CloudType,
+    public CollidingCloudName
 {
 public:

@ -109,8 +112,7 @@ protected:
        // References to the cloud sub-models

            //- Collision model
-            autoPtr<CollisionModel<CollidingCloud<CloudType>>>
-                collisionModel_;
+            autoPtr<CollisionModel<CollidingCloud<CloudType>>> collisionModel_;


        // Initialisation
@ -142,10 +144,10 @@ public:
        CollidingCloud
        (
            const word& cloudName,
+            const dimensionedVector& g,
            const volScalarField& rho,
            const volVectorField& U,
            const volScalarField& mu,
-            const dimensionedVector& g,
            bool readFields = true
        );

@ -184,7 +186,7 @@ public:


    //- Destructor
-    virtual ~CollidingCloud();
+    virtual ~CollidingCloud() = default;


    // Member Functions
--- a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloudName.C
+++ b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloudName.C
@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "CollidingCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(CollidingCloudName, 0);
+}
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@ -30,6 +30,7 @@ License
 #include "integrationScheme.H"
 #include "interpolation.H"
 #include "subCycleTime.H"
+#include "SLGThermo.H"

 #include "InjectionModelList.H"
 #include "DispersionModel.H"
@ -342,15 +343,14 @@ template<class CloudType>
 Foam::KinematicCloud<CloudType>::KinematicCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
    const volScalarField& mu,
-    const dimensionedVector& g,
    bool readFields
 )
 :
-    CloudType(rho.mesh(), cloudName, false),
-    kinematicCloud(),
+    CloudType(cloudName, g, rho, U, mu, false),
    cloudCopyPtr_(nullptr),
    mesh_(rho.mesh()),
    particleProperties_
@ -480,6 +480,29 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
 }


+template<class CloudType>
+Foam::KinematicCloud<CloudType>::KinematicCloud
+(
+    const word& cloudName,
+    const dimensionedVector& g,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    KinematicCloud<CloudType>
+    (
+        cloudName,
+        g,
+        rho,
+        U,
+        thermo.thermo().mu(),
+        readFields
+    )
+{}
+
+
 template<class CloudType>
 Foam::KinematicCloud<CloudType>::KinematicCloud
 (
@ -488,7 +511,6 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
 )
 :
    CloudType(c.mesh_, name, c),
-    kinematicCloud(),
    cloudCopyPtr_(nullptr),
    mesh_(c.mesh_),
    particleProperties_(c.particleProperties_),
@ -561,8 +583,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
    const KinematicCloud<CloudType>& c
 )
 :
-    CloudType(mesh, name, IDLList<parcelType>()),
-    kinematicCloud(),
+    CloudType(mesh, name, c),
    cloudCopyPtr_(nullptr),
    mesh_(mesh),
    particleProperties_
@ -687,8 +708,8 @@ void Foam::KinematicCloud<CloudType>::restoreState()
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::resetSourceTerms()
 {
-    UTrans().field() = Zero;
-    UCoeff().field() = 0.0;
+    UTransRef().field() = Zero;
+    UCoeffRef().field() = 0.0;
 }


@ -725,8 +746,8 @@ void Foam::KinematicCloud<CloudType>::relaxSources
    const KinematicCloud<CloudType>& cloudOldTime
 )
 {
-    this->relax(UTrans_(), cloudOldTime.UTrans(), "U");
-    this->relax(UCoeff_(), cloudOldTime.UCoeff(), "U");
+    this->relax(UTransRef(), cloudOldTime.UTrans_(), "U");
+    this->relax(UCoeffRef(), cloudOldTime.UCoeff_(), "U");
 }


--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@ -59,7 +59,6 @@ SourceFiles

 #include "particle.H"
 #include "Cloud.H"
-#include "kinematicCloud.H"
 #include "IOdictionary.H"
 #include "autoPtr.H"
 #include "Random.H"
@ -80,6 +79,8 @@ namespace Foam

 class integrationScheme;

+class SLGThermo;
+
 template<class CloudType>
 class InjectionModelList;

@ -105,6 +106,8 @@ template<class CloudType>
 class IsotropyModel;


+TemplateName(KinematicCloud);
+
 /*---------------------------------------------------------------------------*\
                       Class KinematicCloud Declaration
 \*---------------------------------------------------------------------------*/
@ -113,7 +116,7 @@ template<class CloudType>
 class KinematicCloud
 :
    public CloudType,
-    public kinematicCloud
+    public KinematicCloudName
 {
 public:

@ -309,15 +312,24 @@ public:


    // Constructors
+        KinematicCloud
+        (
+            const word& cloudName,
+            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const volScalarField& mu,
+            bool readFields = true
+        );

        //- Construct given carrier gas fields
        KinematicCloud
        (
            const word& cloudName,
+            const dimensionedVector& g,
            const volScalarField& rho,
            const volVectorField& U,
-            const volScalarField& mu,
-            const dimensionedVector& g,
+            const SLGThermo& thermo,
            bool readFields = true
        );

@ -527,33 +539,27 @@ public:
                // Momentum

                    //- Return reference to momentum source
-                    inline volVectorField::Internal& UTrans();
+                    inline volVectorField::Internal& UTransRef();

                    //- Return const reference to momentum source
-                    inline const volVectorField::Internal&
-                        UTrans() const;
+                    inline tmp<volVectorField::Internal> UTrans() const;

                     //- Return coefficient for carrier phase U equation
-                    inline volScalarField::Internal& UCoeff();
+                    inline volScalarField::Internal& UCoeffRef();

                    //- Return const coefficient for carrier phase U equation
-                    inline const volScalarField::Internal&
-                        UCoeff() const;
+                    inline tmp<volScalarField::Internal> UCoeff() const;

                    //- Return tmp momentum source term (compressible)
-                    inline tmp<fvVectorMatrix> SU
-                    (
-                        volVectorField& U,
-                        bool incompressible = false
-                    ) const;
+                    inline tmp<fvVectorMatrix> SU(volVectorField& U) const;


        // Check

            //- Total number of parcels
-            virtual label nParcels() const
+            label nParcels() const
            {
-                return CloudType::nParcels();
+                return this->size();
            }

            //- Total mass in system
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@ -436,14 +436,14 @@ Foam::KinematicCloud<CloudType>::cellLengthScale() const

 template<class CloudType>
 inline Foam::DimensionedField<Foam::vector, Foam::volMesh>&
-Foam::KinematicCloud<CloudType>::UTrans()
+Foam::KinematicCloud<CloudType>::UTransRef()
 {
    return *UTrans_;
 }


 template<class CloudType>
-inline const Foam::DimensionedField<Foam::vector, Foam::volMesh>&
+inline Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh>>
 Foam::KinematicCloud<CloudType>::UTrans() const
 {
    return *UTrans_;
@ -452,14 +452,14 @@ Foam::KinematicCloud<CloudType>::UTrans() const

 template<class CloudType>
 inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::KinematicCloud<CloudType>::UCoeff()
+Foam::KinematicCloud<CloudType>::UCoeffRef()
 {
    return *UCoeff_;
 }


 template<class CloudType>
-inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
 Foam::KinematicCloud<CloudType>::UCoeff() const
 {
    return *UCoeff_;
@ -468,7 +468,7 @@ Foam::KinematicCloud<CloudType>::UCoeff() const

 template<class CloudType>
 inline Foam::tmp<Foam::fvVectorMatrix>
-Foam::KinematicCloud<CloudType>::SU(volVectorField& U, bool incompressible)
+Foam::KinematicCloud<CloudType>::SU(volVectorField& U)
 const
 {
    if (debug)
@ -479,38 +479,28 @@ const
            << max(UCoeff()).value() << endl;
    }

-    dimensionSet dim(dimForce);
-    if (incompressible)
-    {
-        dim.reset(dimForce/dimDensity);
-    }
-
    if (solution_.coupled())
    {
+        const auto& deltaT = this->db().time().deltaT();
+
        if (solution_.semiImplicit("U"))
        {
-            volScalarField::Internal
-                Vdt(mesh_.V()*this->db().time().deltaT());
-
-            if (incompressible)
-            {
-                Vdt.dimensions() *= dimDensity;
-            }
+            volScalarField::Internal Vdt(mesh_.V()*deltaT);

            return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;
        }
        else
        {
-            tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dim));
+            auto tfvm = tmp<fvVectorMatrix>::New(U, dimForce);
            fvVectorMatrix& fvm = tfvm.ref();

-            fvm.source() = -UTrans()/(this->db().time().deltaT());
+            fvm.source() = -UTrans()/deltaT;

            return tfvm;
        }
    }

-    return tmp<fvVectorMatrix>::New(U, dim);
+    return tmp<fvVectorMatrix>::New(U, dimForce);
 }


--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudName.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudName.C
@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "KinematicCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(KinematicCloudName, 0);
+}
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
@ -31,6 +31,7 @@ License
 #include "DampingModel.H"
 #include "IsotropyModel.H"
 #include "TimeScaleModel.H"
+#include "SLGThermo.H"

 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //

@ -70,14 +71,14 @@ template<class CloudType>
 Foam::MPPICCloud<CloudType>::MPPICCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
    const volScalarField& mu,
-    const dimensionedVector& g,
    bool readFields
 )
 :
-    CloudType(cloudName, rho, U, mu, g, false),
+    CloudType(cloudName, g, rho, U, mu, false),
    packingModel_(nullptr),
    dampingModel_(nullptr),
    isotropyModel_(nullptr)
@ -102,6 +103,29 @@ Foam::MPPICCloud<CloudType>::MPPICCloud
 }


+template<class CloudType>
+Foam::MPPICCloud<CloudType>::MPPICCloud
+(
+    const word& cloudName,
+    const dimensionedVector& g,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    MPPICCloud<CloudType>
+    (
+        cloudName,
+        g,
+        rho,
+        U,
+        thermo.thermo().mu(),
+        readFields
+    )
+{}
+
+
 template<class CloudType>
 Foam::MPPICCloud<CloudType>::MPPICCloud
 (
@ -131,13 +155,6 @@ Foam::MPPICCloud<CloudType>::MPPICCloud
 {}


-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::MPPICCloud<CloudType>::~MPPICCloud()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

 template<class CloudType>
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H
@ -41,7 +41,6 @@ SourceFiles
 #ifndef MPPICCloud_H
 #define MPPICCloud_H

-#include "particle.H"
 #include "Cloud.H"
 #include "IOdictionary.H"
 #include "autoPtr.H"
@ -55,6 +54,8 @@ namespace Foam

 // Forward declaration of classes

+class SLGThermo;
+
 template<class CloudType>
 class PackingModel;

@ -64,6 +65,8 @@ class DampingModel;
 template<class CloudType>
 class IsotropyModel;

+TemplateName(MPPICCloud);
+
 /*---------------------------------------------------------------------------*\
                       Class MPPICCloud Declaration
 \*---------------------------------------------------------------------------*/
@ -71,7 +74,8 @@ class IsotropyModel;
 template<class CloudType>
 class MPPICCloud
 :
-    public CloudType
+    public CloudType,
+    public MPPICCloudName
 {
 public:

@ -114,12 +118,10 @@ protected:
            autoPtr<PackingModel<MPPICCloud<CloudType>>> packingModel_;

            //- Damping model
-            autoPtr<DampingModel<MPPICCloud<CloudType>>>
-                dampingModel_;
+            autoPtr<DampingModel<MPPICCloud<CloudType>>> dampingModel_;

            //- Exchange model
-            autoPtr<IsotropyModel<MPPICCloud<CloudType>>>
-                isotropyModel_;
+            autoPtr<IsotropyModel<MPPICCloud<CloudType>>> isotropyModel_;


        // Initialisation
@ -136,10 +138,21 @@ public:
        MPPICCloud
        (
            const word& cloudName,
+            const dimensionedVector& g,
            const volScalarField& rho,
            const volVectorField& U,
            const volScalarField& mu,
+            bool readFields = true
+        );
+
+        //- Construct given carrier gas fields
+        MPPICCloud
+        (
+            const word& cloudName,
            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const SLGThermo& thermo,
            bool readFields = true
        );

@ -178,7 +191,7 @@ public:


    //- Destructor
-    virtual ~MPPICCloud();
+    virtual ~MPPICCloud() = default;


    // Member Functions
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloudName.C
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloudName.C
@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MPPICCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(MPPICCloudName, 0);
+}
+
+// ************************************************************************* //
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@ -91,15 +91,14 @@ template<class CloudType>
 Foam::ReactingCloud<CloudType>::ReactingCloud
 (
    const word& cloudName,
+    const dimensionedVector& g,
    const volScalarField& rho,
    const volVectorField& U,
-    const dimensionedVector& g,
    const SLGThermo& thermo,
    bool readFields
 )
 :
-    CloudType(cloudName, rho, U, g, thermo, false),
-    reactingCloud(),
+    CloudType(cloudName, g, rho, U, thermo, false),
    cloudCopyPtr_(nullptr),
    constProps_(this->particleProperties()),
    compositionModel_(nullptr),
@ -155,7 +154,6 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
 )
 :
    CloudType(c, name),
-    reactingCloud(),
    cloudCopyPtr_(nullptr),
    constProps_(c.constProps_),
    compositionModel_(c.compositionModel_->clone()),
@ -195,7 +193,6 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
 )
 :
    CloudType(mesh, name, c),
-    reactingCloud(),
    cloudCopyPtr_(nullptr),
    constProps_(),
    compositionModel_(c.compositionModel_->clone()),
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Kutalmis Bercin	cfc6c89487	WIP: electricPotential: add space charge density hook	2023-02-24 08:34:08 +00:00
Kutalmis Bercin	b93b8affbe	ENH: ChargeCloud: new space charge density function object	2023-02-24 08:34:08 +00:00
Kutalmis Bercin	c63cdfd1d4	ENH: ChargeCloud: add Coulomb force model	2023-02-24 08:34:08 +00:00
Kutalmis Bercin	75cddfe573	ENH: lagrangian: add electric charge cloud	2023-02-24 08:34:08 +00:00
Kutalmis Bercin	089b517ceb	STYLE: electricPotential: simplify read function	2023-02-24 08:34:06 +00:00
Kutalmis Bercin	575ca6cba5	ENH: electricPotential: add fvOptions hook	2023-02-24 08:33:08 +00:00
Kutalmis Bercin	79517b2d39	ENH: electricPotential: replace operandField with getOrReadField	2023-02-21 11:03:55 +00:00
Kutalmis Bercin	55d5f99047	WIP: TUT: tutorial adjustment for backward compatibility	2023-01-31 16:14:43 +00:00
Kutalmis Bercin	d439440453	WIP: TUT: tutorial adjustments for backward compatibility	2023-01-20 16:56:11 +00:00
Andrew Heather	6c66cc58a6	WIP: Updated lagrangian solvers followinf library changes	2023-01-10 13:35:42 +00:00
Andrew Heather	265b3dd2f8	WIP: deprecated lagrangian coal combustion - models added to reactingMultiphase	2023-01-10 13:35:42 +00:00
Andrew Heather	b17a9128a1	ENH: Added new parcelFoam solver based on kinematicParcelFoam	2023-01-10 13:35:42 +00:00
Andrew Heather	31d4148173	STYLE: lagrangian turbulence library - removed 'basic' prefix	2023-01-10 13:35:42 +00:00
Andrew Heather	13b8a3c02a	STYLE: lagrangian intermediate library - removed 'basic' prefix	2023-01-10 13:35:42 +00:00
Andrew Heather	715831610d	WIP: Added new rhoParcelFoam solver based on reactingParcelFoam	2023-01-10 13:35:42 +00:00
Andrew Heather	ab42d1a1e6	WIP: migrated coalCloudList to parcelCloudModelList	2023-01-10 13:35:42 +00:00
Andrew Heather	13cf489f63	WIP: clouds - moved spray cloud and parcels to main intermediate library	2023-01-10 13:35:42 +00:00
Andrew Heather	3277fb931b	WIP: clouds - updated dependent code	2023-01-10 13:35:40 +00:00
Andrew Heather	9e83d31138	WIP: clouds - refactoring	2023-01-10 13:34:47 +00:00
Andrew Heather	878f1bf773	WIP: Clouds - new run-time selection based on openfoam.org changes	2023-01-10 13:34:47 +00:00
Andrew Heather	84348dc3f5	WIP: clouds - added debug - based on openfoam.org changes	2023-01-10 13:34:47 +00:00
Andrew Heather	1f7c038ad7	DUMMY	2023-01-10 13:34:47 +00:00
Andrew Heather	f20361d44a	WIP: Clouds - new run-time selection based on openfoam.org changes	2023-01-10 13:34:33 +00:00
Andrew Heather	580a003d4d	WIP: Removed base classes	2023-01-10 13:34:22 +00:00
mattijs	d74572cae8	CONFIG: set API level to 2212	2022-12-28 09:33:42 +00:00
Andrew Heather	1cd47839ac	CONFIG: Reset version post-release	2022-12-23 13:53:49 +00:00
				`@ -0,0 +1,2 @@`
				`auto parcels = parcelCloudModelList(g, rhoInf, U, muc);`
				`@ -0,0 +1 @@`
				`regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();`