WIP: electricPotential: add space charge density hook

META-INFO/api-info

@@ -1,2 +1,2 @@
 api=2212
-patch=230110
+patch=0

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
 );

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"

applications/solvers/compressible/rhoPimpleFoam/Make/options

@@ -9,7 +9,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/compressible/rhoSimpleFoam/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

@@ -10,7 +10,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -19,7 +19,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 
 EXE_LIBS = \

applications/solvers/incompressible/pimpleFoam/Make/options

@@ -8,7 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

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

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 \

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 \

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

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 \

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 \

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));

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 \

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"

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
 );

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 \

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

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 \

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();
 

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();
 

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 \

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)
     );
 

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
 );
 

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())

applications/solvers/lagrangian/parcelFoam/Make/files

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

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

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);
+    }

applications/solvers/lagrangian/parcelFoam/createClouds.H

@@ -0,0 +1,2 @@
+auto parcels = parcelCloudModelList(g, rhoInf, U, muc);
+

applications/solvers/lagrangian/parcelFoam/createFieldRefs.H

@@ -0,0 +1 @@
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();

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"
+
+

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);
+

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;
+}
+
+
+// ************************************************************************* //

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
 );

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"
 

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)
             {

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 \

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
 );

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)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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;
+}

applications/solvers/lagrangian/rhoParcelFoam/Make/files

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

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

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);
+    }

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);
+}

applications/solvers/lagrangian/rhoParcelFoam/createClouds.H

@@ -0,0 +1,3 @@
+Info<< "\nConstructing cloud" << endl;
+
+auto parcels = parcelCloudModelList(g, rho, U, slgThermo);

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();

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"

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
+        )
+    );
+}
+
+

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)
+);

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);
+    }
+}

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);
+}
+
+// ************************************************************************* //

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;
+}
+
+
+// ************************************************************************* //

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;
+}
+
+// ************************************************************************* //

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;
+}
+
+
+// ************************************************************************* //

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

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
 );

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"

applications/solvers/lagrangian/sprayFoam/Make/options

@@ -35,7 +35,6 @@ EXE_LIBS = \
     -llagrangian \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -llagrangianSpray \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \

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
 );

applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options

@@ -39,7 +39,6 @@ EXE_LIBS = \
     -llagrangian \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -llagrangianSpray \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \

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"

applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options

@@ -34,7 +34,6 @@ EXE_LIBS = \
     -lthermoTools \
     -llagrangian \
     -llagrangianIntermediate \
-    -llagrangianSpray \
     -llagrangianTurbulence \
     -lspecie \
     -lcompressibleTransportModels \

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"

applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options

@@ -40,7 +40,6 @@ EXE_LIBS = \
     -llagrangian \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -llagrangianSpray \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \

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"

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"

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 \

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()
 );

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 \

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();
 

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();
 

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

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 \

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));

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3DKernel.H

@@ -30,8 +30,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_CGALTriangulation3DKernel_H
-#define Foam_CGALTriangulation3DKernel_H
+#ifndef CGALTriangulation3DKernel_H
+#define CGALTriangulation3DKernel_H
 
 // Silence boost bind deprecation warnings (before CGAL-5.2.1)
 #include "CGAL/version.h"
@@ -54,19 +54,9 @@ Description
 //    #include "CGAL/Robust_circumcenter_traits_3.h"
 //    typedef CGAL::Robust_circumcenter_traits_3<baseK>           K;
 
-#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050500000)
-    // Prior to CGAL-5.5
-
     #include "CGAL/Robust_circumcenter_filtered_traits_3.h"
     typedef CGAL::Robust_circumcenter_filtered_traits_3<baseK>  K;
 
-#else
-
-    #include "CGAL/Robust_weighted_circumcenter_filtered_traits_3.h"
-    typedef CGAL::Robust_weighted_circumcenter_filtered_traits_3<baseK> K;
-
-#endif
-
 #else
 
     // Very robust but expensive kernel

applications/utilities/preProcessing/viewFactorsGen/Make/options

@@ -4,6 +4,7 @@ EXE_INC = \
     -Wno-old-style-cast \
     $(COMP_FLAGS) \
     ${CGAL_INC} \
+    -DCGAL_HEADER_ONLY \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
@@ -11,6 +12,7 @@ EXE_INC = \
 
 
 EXE_LIBS = \
+    /* ${CGAL_LIBS} */ \
     -lfiniteVolume \
     -lsurfMesh \
     -lmeshTools \

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.

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.

src/TurbulenceModels/turbulenceModels/laminar/generalizedNewtonian/generalizedNewtonianViscosityModels/Casson/Casson.C

@@ -79,8 +79,8 @@ bool Foam::laminarModels::generalizedNewtonianViscosityModels::Casson::read
 
     coeffs.readEntry("m", m_);
     coeffs.readEntry("tau0", tau0_);
-    coeffs.readEntry("nuMin", nuMin_);
-    coeffs.readEntry("nuMax", nuMax_);
+    coeffs.readEntry("nuMin_", nuMin_);
+    coeffs.readEntry("nuMax_", nuMax_);
 
     return true;
 }

src/finiteArea/Make/files

@@ -55,6 +55,7 @@ $(basicFaPatchFields)/transform/transformFaPatchScalarField.C
 constraintFaPatchFields = $(faPatchFields)/constraint
 $(constraintFaPatchFields)/empty/emptyFaPatchFields.C
 $(constraintFaPatchFields)/processor/processorFaPatchFields.C
+$(constraintFaPatchFields)/processor/processorFaPatchScalarField.C
 $(constraintFaPatchFields)/wedge/wedgeFaPatchFields.C
 $(constraintFaPatchFields)/wedge/wedgeFaPatchScalarField.C
 $(constraintFaPatchFields)/cyclic/cyclicFaPatchFields.C

src/finiteArea/fields/faPatchFields/basic/coupled/coupledFaPatchField.H

@@ -181,6 +181,13 @@ public:
 
         // Coupled interface functionality
 
+            //- Transform given patch component field
+            virtual void transformCoupleField
+            (
+                solveScalarField& f,
+                const direction cmpt
+            ) const = 0;
+
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (

src/finiteArea/fields/faPatchFields/basic/transform/transformFaPatchFieldsFwd.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2016-2017 Wikki Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -23,48 +23,29 @@ License
     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::sprayCloud
-
-Description
-    Virtual abstract base class for templated SprayCloud
-
-SourceFiles
-    sprayCloud.C
-
 \*---------------------------------------------------------------------------*/
 
-#ifndef sprayCloud_H
-#define sprayCloud_H
+#ifndef transformFaPatchFieldsFwd_H
+#define transformFaPatchFieldsFwd_H
 
-#include "typeInfo.H"
+#include "fieldTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-/*---------------------------------------------------------------------------*\
-                       Class sprayCloud Declaration
-\*---------------------------------------------------------------------------*/
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-class sprayCloud
-{
-public:
+template<class Type> class transformFaPatchField;
 
-    //- Runtime type information
-    TypeName("sprayCloud");
+// typedef transformFaPatchField<scalar> transformFaPatchScalarField;
+// typedef transformFaPatchField<vector> transformFaPatchVectorField;
+// typedef transformFaPatchField<tensor> transformFaPatchTensorField;
 
-    // Constructors
-
-        //- Null constructor
-        sprayCloud() = default;
-
-
-    //- Destructor
-    virtual ~sprayCloud() = default;
-};
+// template<class Type> class transformFaPatchField;
 
+makeFaPatchTypeFieldTypedefs(transform)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteArea/fields/faPatchFields/constraint/cyclic/cyclicFaPatchField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -187,7 +187,20 @@ void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
     transformCoupleField(pnf, cmpt);
 
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
+    if (add)
+    {
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
+        }
+    }
+    else
+    {
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+        }
+    }
 }
 
 
@@ -214,11 +227,21 @@ void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
         pnf[facei + sizeby2] = psiInternal[faceCells[facei]];
     }
 
-    // Transform according to the transformation tensors
-    transformCoupleField(pnf);
-
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
+    if (add)
+    {
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
+        }
+    }
+    else
+    {
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+        }
+    }
 }
 
 

src/finiteArea/fields/faPatchFields/constraint/cyclic/cyclicFaPatchField.H

@@ -164,6 +164,16 @@ public:
 
         // Coupled interface functionality
 
+            //- Transform neighbour field
+            virtual void transformCoupleField
+            (
+                solveScalarField& f,
+                const direction cmpt
+            ) const
+            {
+                cyclicLduInterfaceField::transformCoupleField(f, cmpt);
+            }
+
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C

@@ -130,6 +130,13 @@ Foam::processorFaPatchField<Type>::processorFaPatchField
 {}
 
 
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::processorFaPatchField<Type>::~processorFaPatchField()
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -196,8 +203,6 @@ void Foam::processorFaPatchField<Type>::initInterfaceMatrixUpdate
         commsType,
         this->patch().patchInternalField(psiInternal)()
     );
-
-    const_cast<processorFaPatchField<Type>&>(*this).updatedMatrix() = false;
 }
 
 
@@ -214,28 +219,32 @@ void Foam::processorFaPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
-    if (this->updatedMatrix())
-    {
-        return;
-    }
-
-    const labelUList& faceCells = this->patch().edgeFaces();
-
     solveScalarField pnf
     (
-        procPatch_.receive<solveScalar>(commsType, this->size())
+        procPatch_.receive<solveScalar>(commsType, this->size())()
     );
 
-    if (!std::is_arithmetic<Type>::value)
-    {
-        // Transform non-scalar data according to the transformation tensor
-        transformCoupleField(pnf, cmpt);
-    }
+    // Transform according to the transformation tensor
+    transformCoupleField(pnf, cmpt);
 
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 
-    const_cast<processorFaPatchField<Type>&>(*this).updatedMatrix() = true;
+    const labelUList& edgeFaces = this->patch().edgeFaces();
+
+    if (add)
+    {
+        forAll(edgeFaces, elemI)
+        {
+            result[edgeFaces[elemI]] += coeffs[elemI]*pnf[elemI];
+        }
+    }
+    else
+    {
+        forAll(edgeFaces, elemI)
+        {
+            result[edgeFaces[elemI]] -= coeffs[elemI]*pnf[elemI];
+        }
+    }
 }
 
 
@@ -256,8 +265,6 @@ void Foam::processorFaPatchField<Type>::initInterfaceMatrixUpdate
         commsType,
         this->patch().patchInternalField(psiInternal)()
     );
-
-    const_cast<processorFaPatchField<Type>&>(*this).updatedMatrix() = false;
 }
 
 
@@ -273,25 +280,29 @@ void Foam::processorFaPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
-    if (this->updatedMatrix())
-    {
-        return;
-    }
-
-    const labelUList& faceCells = this->patch().edgeFaces();
-
     Field<Type> pnf
     (
-        procPatch_.receive<Type>(commsType, this->size())
+        procPatch_.receive<Type>(commsType, this->size())()
     );
 
-    // Transform according to the transformation tensor
-    transformCoupleField(pnf);
-
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 
-    const_cast<processorFaPatchField<Type>&>(*this).updatedMatrix() = true;
+    const labelUList& edgeFaces = this->patch().edgeFaces();
+
+    if (add)
+    {
+        forAll(edgeFaces, elemI)
+        {
+            result[edgeFaces[elemI]] += coeffs[elemI]*pnf[elemI];
+        }
+    }
+    else
+    {
+        forAll(edgeFaces, elemI)
+        {
+            result[edgeFaces[elemI]] -= coeffs[elemI]*pnf[elemI];
+        }
+    }
 }
 
 

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_processorFaPatchField_H
-#define Foam_processorFaPatchField_H
+#ifndef processorFaPatchField_H
+#define processorFaPatchField_H
 
 #include "coupledFaPatchField.H"
 #include "processorLduInterfaceField.H"
@@ -139,8 +139,9 @@ public:
         }
 
 
-    //- Destructor
-    ~processorFaPatchField() = default;
+    // Destructor
+
+        ~processorFaPatchField();
 
 
     // Member functions
@@ -168,9 +169,18 @@ public:
             //- Return patch-normal gradient
             virtual tmp<Field<Type>> snGrad() const;
 
-
         // Coupled interface functionality
 
+            //- Transform neighbour field
+            virtual void transformCoupleField
+            (
+                solveScalarField& f,
+                const direction cmpt
+            ) const
+            {
+                processorLduInterfaceField::transformCoupleField(f, cmpt);
+            }
+
             //- Initialise neighbour matrix update
             virtual void initInterfaceMatrixUpdate
             (

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchFields.H

@@ -25,10 +25,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_processorFaPatchFields_H
-#define Foam_processorFaPatchFields_H
+#ifndef processorFaPatchFields_H
+#define processorFaPatchFields_H
 
-#include "processorFaPatchField.H"
+#include "processorFaPatchScalarField.H"
 #include "fieldTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchFieldsFwd.H

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016-2017 Wikki 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef processorFaPatchFieldsFwd_H
+#define processorFaPatchFieldsFwd_H
+
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type> class processorFaPatchField;
+
+makeFaPatchTypeFieldTypedefs(processor)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2019 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 "processorFaPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<>
+void Foam::processorFaPatchField<Foam::scalar>::transformCoupleField
+(
+    solveScalarField& f,
+    const direction cmpt
+) const
+{}
+
+
+template<>
+void Foam::processorFaPatchField<Foam::scalar>::initInterfaceMatrixUpdate
+(
+    solveScalarField& result,
+    const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
+    const solveScalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction,
+    const Pstream::commsTypes commsType
+) const
+{
+    procPatch_.send
+    (
+        commsType,
+        patch().patchInternalField(psiInternal)()
+    );
+}
+
+
+template<>
+void Foam::processorFaPatchField<Foam::scalar>::updateInterfaceMatrix
+(
+    solveScalarField& result,
+    const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
+    const solveScalarField&,
+    const scalarField& coeffs,
+    const direction,
+    const Pstream::commsTypes commsType
+) const
+{
+    solveScalarField pnf
+    (
+        procPatch_.receive<solveScalar>(commsType, this->size())()
+    );
+
+    const labelUList& edgeFaces = patch().edgeFaces();
+
+    if (add)
+    {
+        forAll(edgeFaces, facei)
+        {
+            result[edgeFaces[facei]] += coeffs[facei]*pnf[facei];
+        }
+    }
+    else
+    {
+        forAll(edgeFaces, facei)
+        {
+            result[edgeFaces[facei]] -= coeffs[facei]*pnf[facei];
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.H

@@ -5,7 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -23,54 +24,61 @@ License
     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::reactingHeterogeneousCloud
-
-Description
-    Virtual abstract base class for templated ReactingCloud
-
-SourceFiles
-    reactingHeterogeneousCloud.C
-
 \*---------------------------------------------------------------------------*/
 
-#ifndef reactingHeterogeneousCloud_H
-#define reactingHeterogeneousCloud_H
+#ifndef processorFaPatchScalarField_H
+#define processorFaPatchScalarField_H
 
-#include "typeInfo.H"
+#include "processorFaPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-/*---------------------------------------------------------------------------*\
-                 Class reactingHeterogeneousCloud Declaration
-\*---------------------------------------------------------------------------*/
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-class reactingHeterogeneousCloud
-{
-public:
-
-    //- Runtime type information
-    TypeName("reactingHeterogeneousCloud");
+template<>
+void processorFaPatchField<scalar>::transformCoupleField
+(
+    solveScalarField& f,
+    const direction cmpt
+) const;
 
 
-    // Constructors
-
-        //- Null constructor
-        reactingHeterogeneousCloud() = default;
+template<>
+void processorFaPatchField<scalar>::initInterfaceMatrixUpdate
+(
+    solveScalarField& result,
+    const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
+    const solveScalarField&,
+    const scalarField& coeffs,
+    const direction,
+    const Pstream::commsTypes commsType
+) const;
 
 
-    //- Destructor
-    virtual ~reactingHeterogeneousCloud() = default;
-};
+template<>
+void processorFaPatchField<scalar>::updateInterfaceMatrix
+(
+    solveScalarField& result,
+    const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
+    const solveScalarField&,
+    const scalarField& coeffs,
+    const direction,
+    const Pstream::commsTypes commsType
+) const;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/functionObjects/field/expressions/fvExpressionField.C

@@ -389,7 +389,7 @@ bool Foam::functionObjects::fvExpressionField::read(const dictionary& dict)
     }
 
     autowrite_ = dict.getOrDefault("autowrite", false);
-    store_ = dict.getOrDefault("store", true);
+    store_ = dict.getOrDefault("autowrite", true);
 
     // "dimensions" is optional
     dimensions_.clear();

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 \

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)
         {

src/functionObjects/lagrangian/icoUncoupledKinematicCloud/icoUncoupledKinematicCloud.C

@@ -92,10 +92,10 @@ Foam::functionObjects::icoUncoupledKinematicCloud::icoUncoupledKinematicCloud
     kinematicCloud_
     (
         kinematicCloudName_,
+        g_,
         rho_,
         U_,
-        mu_,
-        g_
+        mu_
     )
 {}
 

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

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

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;