BUG: integratedNonUniformTable: correct offsets (fixes #2614)

- with geometryOrder=1, calculate the edge normals from the adjacent
  faces (area-weighted, inverse distance squared) and also
  use that for the Le() calculation.

  Includes the contributions from processor edge neighbours, so it
  should be consistent on both sides.

  This new method (consider as 'beta') contrasts with the current
  standard method that first calculates area-weighted point normals
  and uses the average of them for the edge normal.

  Enable for testing either with a controlDict OptimisationSwitch entry
  "fa:geometryOrder", or on the command-line:

      solverName -opt-switch=fa:geometryOrder=1

.gitlab/issue_templates/bug.md

@@ -49,7 +49,7 @@
 
 <!--
   Providing details of your set-up can help us identify any issues, e.g.
-  OpenFOAM version : v2212|v2206|v2112|v2106|v2012 etc
+  OpenFOAM version : v2206|v2112|v2106|v2012|v2006 etc
   Operating system : ubuntu|openSUSE|centos etc
   Hardware info    : any info that may help?
   Compiler         : gcc|intel|clang etc

META-INFO/api-info

@@ -1,2 +1,2 @@
-api=2212
-patch=0
+api=2206
+patch=220907

README.md

@@ -40,9 +40,9 @@ Violations of the Trademark are monitored, and will be duly prosecuted.
 
 If OpenFOAM has already been compiled on your system, simply source
 the appropriate `etc/bashrc` or `etc/cshrc` file and get started.
-For example, for the OpenFOAM-v2212 version:
+For example, for the OpenFOAM-v2206 version:
 ```
-source /installation/path/OpenFOAM-v2212/etc/bashrc
+source /installation/path/OpenFOAM-v2206/etc/bashrc
 ```
 
 ## Compiling OpenFOAM
@@ -127,8 +127,8 @@ These 3rd-party sources are normally located in a directory parallel
 to the OpenFOAM directory. For example,
 ```
 /path/parent
-|-- OpenFOAM-v2212
-\-- ThirdParty-v2212
+|-- OpenFOAM-v2206
+\-- ThirdParty-v2206
 ```
 There are, however, many cases where this simple convention is inadequate:
 
@@ -136,7 +136,7 @@ There are, however, many cases where this simple convention is inadequate:
   operating system or cluster installation provides it)
 
 * When we have changed the OpenFOAM directory name to some arbitrary
-  directory name, e.g. openfoam-sandbox2212, etc..
+  directory name, e.g. openfoam-sandbox2206, etc..
 
 * When we would like any additional 3rd party software to be located
   inside of the OpenFOAM directory to ensure that the installation is
@@ -156,9 +156,9 @@ when locating the ThirdParty directory with the following precedence:
 2. PREFIX/ThirdParty-VERSION
    * this corresponds to the traditional approach
 3. PREFIX/ThirdParty-vAPI
-   * allows for an updated value of VERSION, *eg*, `v2212-myCustom`,
+   * allows for an updated value of VERSION, *eg*, `v2206-myCustom`,
      without requiring a renamed ThirdParty. The API value would still
-     be `2212` and the original `ThirdParty-v2212/` would be found.
+     be `2206` and the original `ThirdParty-v2206/` would be found.
 4. PREFIX/ThirdParty-API
    * same as the previous example, but using an unadorned API value.
 5. PREFIX/ThirdParty-common

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/options

@@ -3,7 +3,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/overset/lnInclude \
-    -I$(LIB_SRC)/overset/include/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H

@@ -31,25 +31,3 @@
     Info<< "Reading diffusivity DT\n" << endl;
 
     dimensionedScalar DT("DT", dimViscosity, transportProperties);
-
-    bool oversetPatchErrOutput =
-        simple.dict().getOrDefault("oversetPatchErrOutput", false);
-
-    // Dummy phi for oversetPatchErrOutput
-    tmp<surfaceScalarField> tdummyPhi;
-    if (oversetPatchErrOutput)
-    {
-        tdummyPhi = tmp<surfaceScalarField>::New
-        (
-            IOobject
-            (
-                "dummyPhi",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar(dimless, Zero)
-        );
-    }

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/overLaplacianDyMFoam.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,7 +58,6 @@ Description
 #include "fvOptions.H"
 #include "simpleControl.H"
 #include "dynamicFvMesh.H"
-#include "oversetPatchPhiErr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -100,11 +99,6 @@ int main(int argc, char *argv[])
             fvOptions.constrain(TEqn);
             TEqn.solve();
             fvOptions.correct(T);
-
-            if (oversetPatchErrOutput)
-            {
-                oversetPatchPhiErr(TEqn, tdummyPhi.ref());
-            }
         }
 
         #include "write.H"

applications/solvers/basic/potentialFoam/overPotentialFoam/overPotentialFoam.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2022 OpenCFD Ltd
+    Copyright (C) 2017 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -149,6 +149,7 @@ int main(int argc, char *argv[])
 
     mesh.update();
 
+    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     // Since solver contains no time loop it would never execute
     // function objects so do it ourselves

applications/solvers/combustion/fireFoam/createClouds.H

@@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-reactingCloud parcels
+basicReactingCloud parcels
 (
     "reactingCloud1",
-    g,
     rho,
     U,
+    g,
     slgThermo
 );

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -37,7 +37,7 @@ Description
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "reactingCloud.H"
+#include "basicReactingCloud.H"
 #include "surfaceFilmModel.H"
 #include "pyrolysisModelCollection.H"
 #include "radiationModel.H"

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H

@@ -0,0 +1,4 @@
+bool ddtCorr
+(
+    pimple.dict().getOrDefault("ddtCorr", true)
+);

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H

@@ -69,8 +69,6 @@ mesh.setFluxRequired(p.name());
 
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
-// Create bool field with interpolated cells
-#include "createInterpolatedCells.H"
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/overRhoPimpleDyMFoam.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,6 +43,7 @@ Description
 #include "dynamicFvMesh.H"
 #include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
+#include "bound.H"
 #include "pimpleControl.H"
 #include "pressureControl.H"
 #include "CorrectPhi.H"
@@ -88,8 +89,10 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
+        #include "readControls.H"
         #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
@@ -125,6 +128,7 @@ int main(int argc, char *argv[])
         {
             if (pimple.firstIter() || moveMeshOuterCorrectors)
             {
+
                 // Do any mesh changes
                 mesh.update();
 
@@ -133,22 +137,52 @@ int main(int argc, char *argv[])
                     MRF.update();
 
                     #include "setCellMask.H"
-                    #include "setInterpolatedCells.H"
-                    #include "correctRhoPhiFaceMask.H"
+
+                    const surfaceScalarField faceMaskOld
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+                    );
+
+                    // Zero Uf on old faceMask (H-I)
+                    rhoUf() *= faceMaskOld;
+
+                    surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
+
+                    // Update Uf and phi on new C-I faces
+                    rhoUf() += (1-faceMaskOld)*rhoUfint;
+
+                    // Update Uf boundary
+                    forAll(rhoUf().boundaryField(), patchI)
+                    {
+                        rhoUf().boundaryFieldRef()[patchI] =
+                            rhoUfint.boundaryField()[patchI];
+                    }
+
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & rhoUf();
 
                     if (correctPhi)
                     {
-                        // Corrects flux on separated regions
                         #include "correctPhi.H"
                     }
 
+                    // Zero phi on current H-I
+                    const surfaceScalarField faceMask
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask)
+                    );
+
+                    phi *= faceMask;
+                    U   *= cellMask;
+
                      // Make the fluxes relative to the mesh-motion
                     fvc::makeRelative(phi, rho, U);
 
-                    if (checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
+                }
+
+                if (checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
                 }
             }
 

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H

@@ -25,6 +25,17 @@ surfaceScalarField phiHbyA
     fvc::interpolate(rho)*fvc::flux(HbyA)
 );
 
+if (ddtCorr)
+{
+    surfaceScalarField faceMaskOld
+    (
+        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+    );
+
+    phiHbyA +=
+        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf));
+}
+
 fvc::makeRelative(phiHbyA, rho, U);
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
@@ -123,4 +134,8 @@ if (thermo.dpdt())
     }
 }
 
+surfaceScalarField faceMask
+(
+    localMin<scalar>(mesh).interpolate(cellMask)
+);
 phi *= faceMask;

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H

@@ -0,0 +1,9 @@
+#include "readTimeControls.H"
+
+correctPhi = pimple.dict().getOrDefault("correctPhi", false);
+
+checkMeshCourantNo =
+    pimple.dict().getOrDefault("checkMeshCourantNo", false);
+
+
+ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H

@@ -4,7 +4,7 @@
         sqrt
         (
             2*M_PI*sigma*sqr(aMesh.edgeInterpolation::deltaCoeffs())
-            *mag(aMesh.edgeInterpolation::deltaCoeffs())
+            *aMesh.edgeInterpolation::deltaCoeffs()
             /rhol
         )
     ).value()*runTime.deltaT().value();

applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H

@@ -1,6 +1,3 @@
-// Volume-to surface mapping object
-const volSurfaceMapping vsm(aMesh);
-
 volVectorField U
 (
     IOobject
@@ -29,3 +26,6 @@ volScalarField H
     mesh,
     dimensionedScalar(dimLength, Zero)
 );
+
+// Create volume-to surface mapping object
+volSurfaceMapping vsm(aMesh);

applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H

@@ -1,5 +1,5 @@
-    // Volume-to surface mapping object
-    const volSurfaceMapping vsm(aMesh);
+    // Create volume-to surface mapping object
+    volSurfaceMapping vsm(aMesh);
 
     volScalarField Cvf
     (

applications/solvers/finiteArea/surfactantFoam/createVolFields.H

@@ -1,5 +1,5 @@
-    // Volume-to surface mapping object
-    const volSurfaceMapping vsm(aMesh);
+    // Create volume-to surface mapping object
+    volSurfaceMapping vsm(aMesh);
 
     volScalarField Cvf
     (

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H

@@ -124,6 +124,3 @@ dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));
 
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
-
-// Create bool field with interpolated cells
-#include "createInterpolatedCells.H"

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/overBuoyantPimpleDyMFoam.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,7 +50,6 @@ Description
 #include "CorrectPhi.H"
 #include "cellCellStencilObject.H"
 #include "localMin.H"
-#include "oversetAdjustPhi.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -87,6 +86,9 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
+        #include "readTimeControls.H"
+
+        #include "readControls.H"
         #include "readDyMControls.H"
 
         #include "compressibleCourantNo.H"
@@ -126,14 +128,45 @@ int main(int argc, char *argv[])
                     MRF.update();
 
                     #include "setCellMask.H"
-                    #include "setInterpolatedCells.H"
-                    #include "correctRhoPhiFaceMask.H"
+
+                    const surfaceScalarField faceMaskOld
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+                    );
+
+                    // Zero Uf on old faceMask (H-I)
+                    rhoUf() *= faceMaskOld;
+
+                    //fvc::correctRhoUf(rhoUfint, rho, U, phi);
+                    surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
+
+                    // Update Uf and phi on new C-I faces
+                    rhoUf() += (1-faceMaskOld)*rhoUfint;
+
+                    // Update Uf boundary
+                    forAll(rhoUf().boundaryField(), patchI)
+                    {
+                        rhoUf().boundaryFieldRef()[patchI] =
+                            rhoUfint.boundaryField()[patchI];
+                    }
+
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & rhoUf();
 
                     if (correctPhi)
                     {
                         #include "correctPhi.H"
                     }
 
+                    // Zero phi on current H-I
+                    const surfaceScalarField faceMask
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask)
+                    );
+
+                    phi *= faceMask;
+                    U   *= cellMask;
+
                      // Make the fluxes relative to the mesh-motion
                     fvc::makeRelative(phi, rho, U);
                 }

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H

@@ -21,13 +21,16 @@ surfaceScalarField phiHbyA
     fvc::flux(rho*HbyA) + phig
 );
 
-if (adjustFringe)
+if (ddtCorr)
 {
-    fvc::makeRelative(phiHbyA,rho, U);
-    oversetAdjustPhi(phiHbyA, U);
-    fvc::makeAbsolute(phiHbyA,rho, U);
-}
+    surfaceScalarField faceMaskOld
+    (
+        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+    );
 
+    phiHbyA +=
+        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi));
+}
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
@@ -119,4 +122,8 @@ if (thermo.dpdt())
     }
 }
 
+surfaceScalarField faceMask
+(
+    localMin<scalar>(mesh).interpolate(cellMask)
+);
 phi *= faceMask;

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/readControls.H

@@ -0,0 +1,9 @@
+#include "readTimeControls.H"
+
+correctPhi = pimple.dict().getOrDefault("correctPhi", false);
+
+checkMeshCourantNo =
+    pimple.dict().getOrDefault("checkMeshCourantNo", false);
+
+
+ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -113,19 +113,15 @@ int main(int argc, char *argv[])
 
             forAll(fluidRegions, i)
             {
-                fvMesh& mesh = fluidRegions[i];
-
-                #include "readFluidMultiRegionPIMPLEControls.H"
                 #include "setRegionFluidFields.H"
+                #include "readFluidMultiRegionPIMPLEControls.H"
                 #include "solveFluid.H"
             }
 
             forAll(solidRegions, i)
             {
-                fvMesh& mesh = solidRegions[i];
-
-                #include "readSolidMultiRegionPIMPLEControls.H"
                 #include "setRegionSolidFields.H"
+                #include "readSolidMultiRegionPIMPLEControls.H"
                 #include "solveSolid.H"
             }
 
@@ -137,10 +133,8 @@ int main(int argc, char *argv[])
 
                 forAll(fluidRegions, i)
                 {
-                    fvMesh& mesh = fluidRegions[i];
-
-                    #include "readFluidMultiRegionPIMPLEControls.H"
                     #include "setRegionFluidFields.H"
+                    #include "readFluidMultiRegionPIMPLEControls.H"
                     if (!frozenFlow)
                     {
                         Info<< "\nSolving for fluid region "
@@ -172,24 +166,20 @@ int main(int argc, char *argv[])
 
                     forAll(fluidRegions, i)
                     {
-                        fvMesh& mesh = fluidRegions[i];
-
                         Info<< "\nSolving for fluid region "
                             << fluidRegions[i].name() << endl;
-                        #include "readFluidMultiRegionPIMPLEControls.H"
-                        #include "setRegionFluidFields.H"
-                        frozenFlow = true;
-                        #include "solveFluid.H"
+                       #include "setRegionFluidFields.H"
+                       #include "readFluidMultiRegionPIMPLEControls.H"
+                       frozenFlow = true;
+                       #include "solveFluid.H"
                     }
 
                     forAll(solidRegions, i)
                     {
-                        fvMesh& mesh = solidRegions[i];
-
                         Info<< "\nSolving for solid region "
                             << solidRegions[i].name() << endl;
-                        #include "readSolidMultiRegionPIMPLEControls.H"
                         #include "setRegionSolidFields.H"
+                        #include "readSolidMultiRegionPIMPLEControls.H"
                         #include "solveSolid.H"
                     }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -76,21 +76,17 @@ int main(int argc, char *argv[])
 
         forAll(fluidRegions, i)
         {
-            fvMesh& mesh = fluidRegions[i];
-
             Info<< "\nSolving for fluid region "
                 << fluidRegions[i].name() << endl;
-            #include "readFluidMultiRegionSIMPLEControls.H"
             #include "setRegionFluidFields.H"
+            #include "readFluidMultiRegionSIMPLEControls.H"
             #include "solveFluid.H"
         }
 
         forAll(solidRegions, i)
         {
-            fvMesh& mesh = solidRegions[i];
-
-            #include "readSolidMultiRegionSIMPLEControls.H"
             #include "setRegionSolidFields.H"
+            #include "readSolidMultiRegionSIMPLEControls.H"
             #include "solveSolid.H"
         }
 
@@ -103,10 +99,8 @@ int main(int argc, char *argv[])
 
             forAll(fluidRegions, i)
             {
-                fvMesh& mesh = fluidRegions[i];
-
-                #include "readSolidMultiRegionSIMPLEControls.H"
                 #include "setRegionFluidFields.H"
+                #include "readSolidMultiRegionSIMPLEControls.H"
                 if (!frozenFlow)
                 {
                     #include "pEqn.H"
@@ -127,24 +121,20 @@ int main(int argc, char *argv[])
 
                 forAll(fluidRegions, i)
                 {
-                    fvMesh& mesh = fluidRegions[i];
-
                     Info<< "\nSolving for fluid region "
                         << fluidRegions[i].name() << endl;
-                    #include "readFluidMultiRegionSIMPLEControls.H"
-                    #include "setRegionFluidFields.H"
-                    frozenFlow = true;
-                    #include "solveFluid.H"
+                   #include "setRegionFluidFields.H"
+                   #include "readFluidMultiRegionSIMPLEControls.H"
+                   frozenFlow = true;
+                   #include "solveFluid.H"
                 }
 
                 forAll(solidRegions, i)
                 {
-                    fvMesh& mesh = solidRegions[i];
-
                     Info<< "\nSolving for solid region "
                         << solidRegions[i].name() << endl;
-                    #include "readSolidMultiRegionSIMPLEControls.H"
                     #include "setRegionSolidFields.H"
+                    #include "readSolidMultiRegionSIMPLEControls.H"
                     #include "solveSolid.H"
                 }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H

@@ -5,5 +5,3 @@
 
     const bool momentumPredictor =
         simple.getOrDefault("momentumPredictor", true);
-
-    simple.readIfPresent("frozenFlow", frozenFlowFluid[i]);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -1,3 +1,5 @@
+    const fvMesh& mesh = fluidRegions[i];
+
     rhoThermo& thermo = thermoFluid[i];
     thermo.validate(args.executable(), "h", "e");
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/chtMultiRegionTwoPhaseEulerFoam.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -108,23 +108,19 @@ int main(int argc, char *argv[])
 
             forAll(fluidRegions, i)
             {
-                fvMesh& mesh = fluidRegions[i];
-
                 Info<< "\nSolving for fluid region "
                     << fluidRegions[i].name() << endl;
-                #include "readFluidMultiRegionPIMPLEControls.H"
                 #include "setRegionFluidFields.H"
+                #include "readFluidMultiRegionPIMPLEControls.H"
                 #include "solveFluid.H"
             }
 
             forAll(solidRegions, i)
             {
-                fvMesh& mesh = solidRegions[i];
-
                 Info<< "\nSolving for solid region "
                     << solidRegions[i].name() << endl;
-                #include "readSolidMultiRegionPIMPLEControls.H"
                 #include "setRegionSolidFields.H"
+                #include "readSolidMultiRegionPIMPLEControls.H"
                 #include "solveSolid.H"
             }
 
@@ -139,24 +135,20 @@ int main(int argc, char *argv[])
 
                     forAll(fluidRegions, i)
                     {
-                        fvMesh& mesh = fluidRegions[i];
-
                         Info<< "\nSolving for fluid region "
                             << fluidRegions[i].name() << endl;
-                        #include "readFluidMultiRegionPIMPLEControls.H"
-                        #include "setRegionFluidFields.H"
-                        frozenFlow = true;
-                        #include "solveFluid.H"
+                       #include "setRegionFluidFields.H"
+                       #include "readFluidMultiRegionPIMPLEControls.H"
+                       frozenFlow = true;
+                       #include "solveFluid.H"
                     }
 
                     forAll(solidRegions, i)
                     {
-                        fvMesh& mesh = solidRegions[i];
-
                         Info<< "\nSolving for solid region "
                             << solidRegions[i].name() << endl;
-                        #include "readSolidMultiRegionPIMPLEControls.H"
                         #include "setRegionSolidFields.H"
+                        #include "readSolidMultiRegionPIMPLEControls.H"
                         #include "solveSolid.H"
                     }
                 }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H

@@ -9,5 +9,3 @@
     (
         pimpleDict.getOrDefault<int>("nEnergyCorrectors", 1)
     );
-
-    pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H

@@ -1,3 +1,5 @@
+    fvMesh& mesh = fluidRegions[i];
+
     twoPhaseSystem& fluid  = phaseSystemFluid[i];
 
     phaseModel& phase1 = fluid.phase1();

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H

@@ -8,5 +8,3 @@
 
     const bool momentumPredictor =
         pimple.getOrDefault("momentumPredictor", true);
-
-    pimple.readIfPresent("frozenFlow", frozenFlowFluid[i]);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -1,3 +1,5 @@
+    fvMesh& mesh = fluidRegions[i];
+
     CombustionModel<rhoReactionThermo>& reaction = reactionFluid[i];
 
     rhoReactionThermo& thermo = reaction.thermo();

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -35,7 +35,7 @@
                 (
                     solidRegions[i],
                     thermos[i],
-                    coordinateSystem::typeName
+                    coordinateSystem::typeName_()
                 )
             );
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -1,3 +1,4 @@
+fvMesh& mesh = solidRegions[i];
 solidThermo& thermo = thermos[i];
 
 tmp<volScalarField> trho = thermo.rho();

applications/solvers/heatTransfer/solidFoam/createFields.H

@@ -15,7 +15,7 @@ if (!thermo.isotropic())
     (
         mesh,
         thermo,
-        coordinateSystem::typeName
+        coordinateSystem::typeName_()
     );
 
     tmp<volVectorField> tkappaByCp = thermo.Kappa()/thermo.Cp();

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,10 +37,7 @@ scalar meanCoNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField phiMask
-    (
-        localMin<scalar>(mesh).interpolate(cellMask + interpolatedCells)
-    );
+    surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     scalarField sumPhi(fvc::surfaceSum(mag(phiMask*phi))().internalField());
 

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H

@@ -1,4 +1,5 @@
 // Solve the Momentum equation
+
 MRF.correctBoundaryVelocity(U);
 
 tmp<fvVectorMatrix> tUEqn

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H

@@ -0,0 +1,26 @@
+#include "createTimeControls.H"
+
+bool correctPhi
+(
+    pimple.dict().getOrDefault("correctPhi", false)
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().getOrDefault("checkMeshCourantNo", false)
+);
+
+bool massFluxInterpolation
+(
+    pimple.dict().getOrDefault("massFluxInterpolation", false)
+);
+
+bool adjustFringe
+(
+    pimple.dict().getOrDefault("oversetAdjustPhi", false)
+);
+
+bool ddtCorr
+(
+    pimple.dict().getOrDefault("ddtCorr", true)
+);

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H

@@ -0,0 +1,273 @@
+// Interpolation used
+interpolationCellPoint<vector> UInterpolator(HbyA);
+
+// Determine faces on outside of interpolated cells
+bitSet isOwnerInterpolatedFace(mesh.nInternalFaces());
+bitSet isNeiInterpolatedFace(mesh.nInternalFaces());
+
+// Determine donor cells
+labelListList donorCell(mesh.nInternalFaces());
+
+scalarListList weightCellCells(mesh.nInternalFaces());
+
+// Interpolated HbyA faces
+vectorField UIntFaces(mesh.nInternalFaces(), Zero);
+
+// Determine receptor neighbour cells
+labelList receptorNeigCell(mesh.nInternalFaces(), -1);
+
+{
+    const cellCellStencilObject& overlap = Stencil::New(mesh);
+    const labelList& cellTypes = overlap.cellTypes();
+    const labelIOList& zoneID = overlap.zoneID();
+
+    label nZones = gMax(zoneID)+1;
+    PtrList<fvMeshSubset> meshParts(nZones);
+    labelList nCellsPerZone(nZones, Zero);
+
+    // A mesh subset for each zone
+    forAll(meshParts, zonei)
+    {
+        meshParts.set
+        (
+            zonei,
+            // Select cells where the zoneID == zonei
+            new fvMeshSubset(mesh, zonei, zoneID)
+        );
+    }
+
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        label ownType = cellTypes[mesh.faceOwner()[faceI]];
+        label neiType = cellTypes[mesh.faceNeighbour()[faceI]];
+        if
+        (
+            ownType == cellCellStencil::INTERPOLATED
+            && neiType == cellCellStencil::CALCULATED
+        )
+        {
+            isOwnerInterpolatedFace.set(faceI);
+
+            const vector& fc = mesh.faceCentres()[faceI];
+
+            for (label zoneI = 0; zoneI < nZones; zoneI++)
+            {
+                if (zoneI != zoneID[mesh.faceOwner()[faceI]])
+                {
+                    const fvMesh& partMesh = meshParts[zoneI].subMesh();
+                    const labelList& cellMap = meshParts[zoneI].cellMap();
+                    label cellI = partMesh.findCell(fc);
+
+                    if (cellI != -1)
+                    {
+                         // Determine weights
+                        labelList stencil(partMesh.cellCells()[cellI]);
+
+                        stencil.append(cellI);
+
+                        label st = stencil.size();
+
+                        donorCell[faceI].setSize(st);
+
+                        weightCellCells[faceI].setSize(st);
+
+                        scalarField weights(st);
+
+                        forAll(stencil, i)
+                        {
+                            scalar d = mag
+                            (
+                                partMesh.cellCentres()[stencil[i]]
+                              - fc
+                            );
+                            weights[i] = 1.0/d;
+                            donorCell[faceI][i] = cellMap[stencil[i]];
+                        }
+                        weights /= sum(weights);
+
+                        weightCellCells[faceI] = weights;
+
+                        forAll(stencil, i)
+                        {
+                            UIntFaces[faceI] +=
+                                weightCellCells[faceI][i]
+                               *UInterpolator.interpolate
+                                (
+                                    fc,
+                                    donorCell[faceI][i]
+                                );
+                        }
+
+                        break;
+                    }
+                }
+            }
+
+            receptorNeigCell[faceI] = mesh.faceNeighbour()[faceI];
+        }
+        else if
+        (
+            ownType == cellCellStencil::CALCULATED
+            && neiType == cellCellStencil::INTERPOLATED
+        )
+        {
+            isNeiInterpolatedFace.set(faceI);
+
+            const vector& fc = mesh.faceCentres()[faceI];
+            for (label zoneI = 0; zoneI < nZones; zoneI++)
+            {
+                if (zoneI != zoneID[mesh.faceNeighbour()[faceI]])
+                {
+                    const fvMesh& partMesh = meshParts[zoneI].subMesh();
+                    const labelList& cellMap = meshParts[zoneI].cellMap();
+                    label cellI = partMesh.findCell(fc);
+
+                    if (cellI != -1)
+                    {
+                        // Determine weights
+                        labelList stencil(partMesh.cellCells()[cellI]);
+
+                        stencil.append(cellI);
+
+                        label st = stencil.size();
+
+                        donorCell[faceI].setSize(st);
+
+                        weightCellCells[faceI].setSize(st);
+
+                        scalarField weights(st);
+
+                        forAll(stencil, i)
+                        {
+                            scalar d = mag
+                            (
+                                partMesh.cellCentres()[stencil[i]]
+                              - fc
+                            );
+                            weights[i] = 1.0/d;
+                            donorCell[faceI][i] = cellMap[stencil[i]];
+                        }
+                        weights /= sum(weights);
+
+                        weightCellCells[faceI] = weights;
+
+                        forAll(stencil, i)
+                        {
+                            UIntFaces[faceI] +=
+                                weightCellCells[faceI][i]
+                               *UInterpolator.interpolate
+                                (
+                                    fc,
+                                    donorCell[faceI][i]
+                                );
+                        }
+
+                        break;
+                    }
+                }
+            }
+
+            receptorNeigCell[faceI] = mesh.faceOwner()[faceI];
+        }
+    }
+}
+
+// contravariant U
+vectorField U1Contrav(mesh.nInternalFaces(), Zero);
+
+surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf());
+
+forAll(isNeiInterpolatedFace, faceI)
+{
+    label cellId = -1;
+    if (isNeiInterpolatedFace.test(faceI))
+    {
+        cellId = mesh.faceNeighbour()[faceI];
+    }
+    else if (isOwnerInterpolatedFace.test(faceI))
+    {
+        cellId = mesh.faceOwner()[faceI];
+    }
+
+    if (cellId != -1)
+    {
+        const vector& n = faceNormals[faceI];
+        vector n1(Zero);
+
+        // 2-D cases
+        if (mesh.nSolutionD() == 2)
+        {
+            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+            {
+                if (mesh.geometricD()[cmpt] == -1)
+                {
+                    switch (cmpt)
+                    {
+                        case vector::X:
+                        {
+                            n1 = vector(0, n.z(), -n.y());
+                            break;
+                        }
+
+                        case vector::Y:
+                        {
+                            n1 = vector(n.z(), 0, -n.x());
+                            break;
+                        }
+
+                        case vector::Z:
+                        {
+                            n1 = vector(n.y(), -n.x(), 0);
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+        else if (mesh.nSolutionD() == 3)
+        {
+            //Determine which is the primary direction
+            if (mag(n.x()) > mag(n.y()) && mag(n.x()) > mag(n.z()))
+            {
+                n1 = vector(n.y(), -n.x(), 0);
+            }
+            else if (mag(n.y()) > mag(n.z()))
+            {
+                n1 = vector(0, n.z(), -n.y());
+            }
+            else
+            {
+                n1 = vector(-n.z(), 0, n.x());
+            }
+        }
+        n1.normalise();
+
+        const vector n2 = normalised(n ^ n1);
+
+        tensor rot =
+            tensor
+            (
+               n.x() ,n.y(), n.z(),
+               n1.x() ,n1.y(), n1.z(),
+               n2.x() ,n2.y(), n2.z()
+            );
+
+//         tensor rot =
+//             tensor
+//             (
+//                n  & x ,n  & y, n  & z,
+//                n1 & x ,n1 & y, n1 & z,
+//                n2 & x ,n2 & y, n2 & z
+//             );
+
+        U1Contrav[faceI].x() =
+           2*transform(rot, UIntFaces[faceI]).x()
+           - transform(rot, HbyA[receptorNeigCell[faceI]]).x();
+
+        U1Contrav[faceI].y() = transform(rot, HbyA[cellId]).y();
+
+        U1Contrav[faceI].z() = transform(rot, HbyA[cellId]).z();
+
+        HbyA[cellId] = transform(inv(rot), U1Contrav[faceI]);
+    }
+}

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/overPimpleDyMFoam.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,9 +46,12 @@ Description
 #include "fvOptions.H"
 
 #include "cellCellStencilObject.H"
+#include "zeroGradientFvPatchFields.H"
 #include "localMin.H"
+#include "interpolationCellPoint.H"
+#include "transform.H"
+#include "fvMeshSubset.H"
 #include "oversetAdjustPhi.H"
-#include "oversetPatchPhiErr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,9 +68,10 @@ int main(int argc, char *argv[])
     #include "setRootCaseLists.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "createDyMControls.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     #include "createFields.H"
     #include "createUf.H"
     #include "createMRF.H"
@@ -84,9 +88,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readDyMControls.H"
-        #include "readOversetDyMControls.H"
-
+        #include "readControls.H"
         #include "CourantNo.H"
 
         #include "setDeltaT.H"
@@ -95,20 +97,45 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        mesh.update();
+        bool changed = mesh.update();
 
-        if (mesh.changing())
+        if (changed)
         {
             #include "setCellMask.H"
             #include "setInterpolatedCells.H"
-            #include "correctPhiFaceMask.H"
 
-            fvc::makeRelative(phi, U);
+            surfaceScalarField faceMaskOld
+            (
+                localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+            );
 
-            if (checkMeshCourantNo)
-            {
-                #include "meshCourantNo.H"
-            }
+            // Zero Uf on old faceMask (H-I)
+            Uf *= faceMaskOld;
+            // Update Uf and phi on new C-I faces
+            Uf += (1-faceMaskOld)*fvc::interpolate(U);
+            phi = mesh.Sf() & Uf;
+
+            // Zero phi on current H-I
+            surfaceScalarField faceMask
+            (
+                localMin<scalar>(mesh).interpolate(cellMask)
+            );
+            phi *= faceMask;
+        }
+
+
+        if (mesh.changing() && correctPhi)
+        {
+            // Calculate absolute flux from the mapped surface velocity
+            #include "correctPhi.H"
+        }
+
+        // Make the flux relative to the mesh motion
+        fvc::makeRelative(phi, U);
+
+        if (mesh.changing() && checkMeshCourantNo)
+        {
+            #include "meshCourantNo.H"
         }
 
         // --- Pressure-velocity PIMPLE corrector loop

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H

@@ -1,11 +1,36 @@
+// Option 1: interpolate rAU, do not block out rAU on blocked cells
 volScalarField rAU("rAU", 1.0/UEqn.A());
+mesh.interpolate(rAU);
 
+// Option 2: do not interpolate rAU but block out rAU
+//surfaceScalarField rAUf("rAUf", fvc::interpolate(blockedCells*rAU));
+
+
+// Option 3: do not interpolate rAU but zero out rAUf on faces on holes
+// But what about:
+//
+//   H
+// H I C C C C
+//   H
+//
 surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 volVectorField H("H", UEqn.H());
 
 volVectorField HbyA("HbyA", U);
 HbyA = constrainHbyA(rAU*H, U, p);
 
+if (massFluxInterpolation)
+{
+    #include "interpolatedFaces.H"
+}
+
+if (runTime.outputTime())
+{
+    H.write();
+    rAU.write();
+    HbyA.write();
+}
+
 if (pimple.nCorrPISO() <= 1)
 {
     tUEqn.clear();
@@ -13,16 +38,33 @@ if (pimple.nCorrPISO() <= 1)
 
 phiHbyA = fvc::flux(HbyA);
 
+if (ddtCorr)
+{
+    surfaceScalarField faceMaskOld
+    (
+        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+    );
+    phiHbyA += rAUf*faceMaskOld*fvc::ddtCorr(U, Uf);
+}
+
 MRF.makeRelative(phiHbyA);
 
+// WIP
+if (p.needReference())
+{
+    fvc::makeRelative(phiHbyA, U);
+    adjustPhi(phiHbyA, U, p);
+    fvc::makeAbsolute(phiHbyA, U);
+}
+
+
 if (adjustFringe)
 {
     fvc::makeRelative(phiHbyA, U);
-    oversetAdjustPhi(phiHbyA, U, zoneIdMass);
+    oversetAdjustPhi(phiHbyA, U);
     fvc::makeAbsolute(phiHbyA, U);
 }
 
-
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
@@ -37,26 +79,27 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA - pEqn.flux();
-        pEqn.relax();
-        U =
-            cellMask*
-            (
-                HbyA
-              - rAU*fvc::reconstruct((pEqn.flux())/rAUf)
-            );
-        U.correctBoundaryConditions();
-        fvOptions.correct(U);
-    }
-
-    if (oversetPatchErrOutput)
-    {
-        oversetPatchPhiErr(pEqn, phiHbyA);
+        // option 2:
+        // rAUf*fvc::snGrad(p)*mesh.magSf();
     }
 }
 
-// Excludes error in interpolated/hole cells
+
 #include "continuityErrs.H"
 
+// Explicitly relax pressure for momentum corrector
+p.relax();
+volVectorField gradP(fvc::grad(p));
+
+// Option 2: zero out velocity on blocked out cells
+//U = HbyA - rAU*cellMask*gradP;
+// Option 3: zero out velocity on blocked out cells
+// This is needed for the scalar Eq (k,epsilon, etc)
+// which can use U as source term
+U = cellMask*(HbyA - rAU*gradP);
+U.correctBoundaryConditions();
+
+fvOptions.correct(U);
 
 {
     Uf = fvc::interpolate(U);
@@ -66,4 +109,9 @@ while (pimple.correctNonOrthogonal())
 
 // Make the fluxes relative to the mesh motion
 fvc::makeRelative(phi, U);
+
+surfaceScalarField faceMask
+(
+    localMin<scalar>(mesh).interpolate(cellMask)
+);
 phi *= faceMask;

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/readControls.H

@@ -0,0 +1,10 @@
+#include "readTimeControls.H"
+
+correctPhi = pimple.dict().getOrDefault("correctPhi", false);
+
+checkMeshCourantNo = pimple.dict().getOrDefault("checkMeshCourantNo", false);
+
+massFluxInterpolation =
+    pimple.dict().getOrDefault("massFluxInterpolation", false);
+
+ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

applications/solvers/incompressible/simpleFoam/overSimpleFoam/createOversetFields.H

@@ -24,3 +24,7 @@ bool adjustFringe
 (
     simple.dict().getOrDefault("oversetAdjustPhi", false)
 );
+bool massFluxInterpolation
+(
+    simple.dict().getOrDefault("massFluxInterpolation", false)
+);

applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H

@@ -1,10 +1,18 @@
 {
+    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
 
+    //mesh.interpolate(HbyA);
+    if (massFluxInterpolation)
+    {
+        #include "interpolatedFaces.H"
+    }
+
     tUEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C

@@ -42,11 +42,11 @@ Description
 #include "CorrectPhi.H"
 
 #ifdef MPPIC
-    #include "kinematicCloud.H"
-    #define kinematicTypeCloud kinematicCloud
+    #include "basicKinematicCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCloud
 #else
-    #include "kinematicCollidingCloud.H"
-    #define kinematicTypeCloud kinematicCollidingCloud
+    #include "basicKinematicCollidingCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCollidingCloud
 #endif
 
 int main(int argc, char *argv[])
@@ -88,7 +88,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // Store the particle positions
-        kCloud.storeGlobalPositions();
+        kinematicCloud.storeGlobalPositions();
 
         mesh.update();
 
@@ -111,15 +111,15 @@ int main(int argc, char *argv[])
         continuousPhaseTransport.correct();
         muc = rhoc*continuousPhaseTransport.nu();
 
-        kCloud.evolve();
+        kinematicCloud.evolve();
 
         // Update continuous phase volume fraction field
-        alphac = max(1.0 - kCloud.theta(), alphacMin);
+        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
         alphac.correctBoundaryConditions();
         alphacf = fvc::interpolate(alphac);
         alphaPhic = alphacf*phic;
 
-        fvVectorMatrix cloudSU(kCloud.SU(Uc));
+        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
         volVectorField cloudVolSUSu
         (
             IOobject

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options

@@ -12,8 +12,6 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options

@@ -11,11 +11,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \

applications/solvers/lagrangian/DPMFoam/DPMFoam.C

@@ -43,11 +43,11 @@ Description
 #include "pimpleControl.H"
 
 #ifdef MPPIC
-    #include "kinematicCloud.H"
-    #define kinematicTypeCloud kinematicCloud
+    #include "basicKinematicCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCloud
 #else
-    #include "kinematicCollidingCloud.H"
-    #define kinematicTypeCloud kinematicCollidingCloud
+    #include "basicKinematicCollidingCloud.H"
+    #define basicKinematicTypeCloud basicKinematicCollidingCloud
 #endif
 
 int main(int argc, char *argv[])
@@ -91,16 +91,16 @@ int main(int argc, char *argv[])
         continuousPhaseTransport.correct();
         muc = rhoc*continuousPhaseTransport.nu();
 
-        Info<< "Evolving " << kCloud.name() << endl;
-        kCloud.evolve();
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+        kinematicCloud.evolve();
 
         // Update continuous phase volume fraction field
-        alphac = max(1.0 - kCloud.theta(), alphacMin);
+        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
         alphac.correctBoundaryConditions();
         alphacf = fvc::interpolate(alphac);
         alphaPhic = alphacf*phic;
 
-        fvVectorMatrix cloudSU(kCloud.SU(Uc));
+        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
         volVectorField cloudVolSUSu
         (
             IOobject

applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options

@@ -10,8 +10,6 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \

applications/solvers/lagrangian/DPMFoam/Make/options

@@ -10,11 +10,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \

applications/solvers/lagrangian/DPMFoam/createFields.H

@@ -125,13 +125,13 @@ const word kinematicCloudName
 );
 
 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-kinematicTypeCloud kCloud
+basicKinematicTypeCloud kinematicCloud
 (
     kinematicCloudName,
-    g,
     rhoc,
     Uc,
-    muc
+    muc,
+    g
 );
 
 // Particle fraction upper limit
@@ -139,13 +139,13 @@ scalar alphacMin
 (
     1.0
   - (
-        kCloud.particleProperties().subDict("constantProperties")
+        kinematicCloud.particleProperties().subDict("constantProperties")
        .get<scalar>("alphaMax")
     )
 );
 
 // Update alphac from the particle locations
-alphac = max(1.0 - kCloud.theta(), alphacMin);
+alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
 alphac.correctBoundaryConditions();
 
 surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));

applications/solvers/lagrangian/coalChemistryFoam/Make/options

@@ -6,6 +6,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@@ -34,6 +35,7 @@ EXE_LIBS = \
     -llagrangian \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
+    -lcoalCombustion\
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -37,8 +37,8 @@ Description
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "thermoCloud.H"
-#include "reactingMultiphaseCloud.H"
+#include "basicThermoCloud.H"
+#include "coalCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "fvOptions.H"

applications/solvers/lagrangian/coalChemistryFoam/createClouds.H

@@ -1,19 +1,19 @@
 Info<< "\nConstructing coal cloud" << endl;
-reactingMultiphaseCloud coalParcels
+coalCloud coalParcels
 (
     "coalCloud1",
-    g,
     rho,
     U,
+    g,
     slgThermo
 );
 
 Info<< "\nConstructing limestone cloud" << endl;
-thermoCloud limestoneParcels
+basicThermoCloud limestoneParcels
 (
     "limestoneCloud1",
-    g,
     rho,
     U,
+    g,
     slgThermo
 );

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options

@@ -9,8 +9,6 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/createFields.H

@@ -63,13 +63,13 @@ const word kinematicCloudName
 );
 
 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-kinematicCollidingCloud kCloud
+basicKinematicCollidingCloud kinematicCloud
 (
     kinematicCloudName,
-    g,
     rhoInf,
     U,
-    mu
+    mu,
+    g
 );
 
 IOobject Hheader

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options

@@ -10,8 +10,6 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C

@@ -41,7 +41,7 @@ Description
 #include "dynamicFvMesh.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
-#include "kinematicCollidingCloud.H"
+#include "basicKinematicCollidingCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -76,19 +76,19 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        kCloud.storeGlobalPositions();
+        kinematicCloud.storeGlobalPositions();
 
         mesh.update();
 
         U.correctBoundaryConditions();
 
-        Info<< "Evolving " << kCloud.name() << endl;
+        Info<< "Evolving " << kinematicCloud.name() << endl;
 
         laminarTransport.correct();
 
         mu = laminarTransport.nu()*rhoInfValue;
 
-        kCloud.evolve();
+        kinematicCloud.evolve();
 
         runTime.write();
 

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelFoam.C

@@ -40,7 +40,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
-#include "kinematicCollidingCloud.H"
+#include "basicKinematicCollidingCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,13 +75,13 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        Info<< "Evolving " << kCloud.name() << endl;
+        Info<< "Evolving " << kinematicCloud.name() << endl;
 
         laminarTransport.correct();
 
         mu = laminarTransport.nu()*rhoInfValue;
 
-        kCloud.evolve();
+        kinematicCloud.evolve();
 
         runTime.write();
 

applications/solvers/lagrangian/kinematicParcelFoam/Make/options

@@ -8,12 +8,6 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \

applications/solvers/lagrangian/kinematicParcelFoam/UEqn.H

@@ -6,7 +6,7 @@
       + MRF.DDt(U)
       + turbulence->divDevReff(U)
      ==
-        scalar(1)/rhoInfValue*parcels.SU(U)
+        parcels.SU(U, true)
       + fvOptions(U)
     );
 

applications/solvers/lagrangian/kinematicParcelFoam/createClouds.H

@@ -4,12 +4,13 @@ const word kinematicCloudName
 );
 
 Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
-kinematicCloud parcels
+
+basicKinematicCloud parcels
 (
     kinematicCloudName,
-    g,
     rhoInf,
     U,
-    muc
+    muc,
+    g
 );
 

applications/solvers/lagrangian/kinematicParcelFoam/kinematicParcelFoam.C

@@ -40,7 +40,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "surfaceFilmModel.H"
-#include "kinematicCloud.H"
+#include "basicKinematicCloud.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
 #include "CorrectPhi.H"
@@ -67,6 +67,7 @@ int main(int argc, char *argv[])
     #include "createDyMControls.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
+    #include "createRegionControls.H"
     #include "createUfIfPresent.H"
 
     turbulence->validate();
@@ -94,7 +95,7 @@ int main(int argc, char *argv[])
         // Do any mesh changes
         mesh.update();
 
-        if (pimple.solveFlow() && mesh.changing())
+        if (solvePrimaryRegion && mesh.changing())
         {
             MRF.update();
 
@@ -119,7 +120,7 @@ int main(int argc, char *argv[])
         parcels.evolve();
         surfaceFilm.evolve();
 
-        if (pimple.solveFlow())
+        if (solvePrimaryRegion)
         {
             // --- PIMPLE loop
             while (pimple.loop())

applications/solvers/lagrangian/parcelFoam/Make/files

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

applications/solvers/lagrangian/parcelFoam/UEqn.H

@@ -1,22 +0,0 @@
-    MRF.correctBoundaryVelocity(U);
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(U) + fvm::div(phi, U)
-      + MRF.DDt(U)
-      + turbulence->divDevReff(U)
-     ==
-        scalar(1)/rhoInfValue*parcels.SU(U)
-      + fvOptions(U)
-    );
-
-    UEqn.relax();
-
-    fvOptions.constrain(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve(UEqn == -fvc::grad(p));
-
-        fvOptions.correct(U);
-    }

applications/solvers/lagrangian/parcelFoam/createClouds.H

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

applications/solvers/lagrangian/parcelFoam/createFieldRefs.H

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

applications/solvers/lagrangian/parcelFoam/createFields.H

@@ -1,85 +0,0 @@
-#include "readGravitationalAcceleration.H"
-
-Info<< "Reading field p\n" << endl;
-volScalarField p
-(
-    IOobject
-    (
-        "p",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-Info<< "\nReading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "createPhi.H"
-
-singlePhaseTransportModel laminarTransport(U, phi);
-
-dimensionedScalar rhoInfValue
-(
-    "rhoInf",
-    dimDensity,
-    laminarTransport
-);
-
-volScalarField rhoInf
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    mesh,
-    rhoInfValue
-);
-
-volScalarField muc
-(
-    IOobject
-    (
-        "muc",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    rhoInf*laminarTransport.nu()
-);
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<incompressible::turbulenceModel> turbulence
-(
-    incompressible::turbulenceModel::New(U, phi, laminarTransport)
-);
-
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, pimple.dict(), pRefCell, pRefValue);
-mesh.setFluxRequired(p.name());
-
-#include "createMRF.H"
-#include "createClouds.H"
-#include "createSurfaceFilmModel.H"
-#include "createFvOptions.H"
-
-

applications/solvers/lagrangian/parcelFoam/pEqn.H

@@ -1,58 +0,0 @@
-
-volScalarField rAU(1.0/UEqn.A());
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
-
-if (pimple.ddtCorr())
-{
-    phiHbyA += MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf));
-}
-
-MRF.makeRelative(phiHbyA);
-
-if (p.needReference())
-{
-    fvc::makeRelative(phiHbyA, U);
-    adjustPhi(phiHbyA, U, p);
-    fvc::makeAbsolute(phiHbyA, U);
-}
-
-
-// Update the pressure BCs to ensure flux consistency
-constrainPressure(p, U, phiHbyA, rAU, MRF);
-
-// Non-orthogonal pressure corrector loop
-while (pimple.correctNonOrthogonal())
-{
-    fvScalarMatrix pEqn
-    (
-        fvm::laplacian(rAU, p)
-     ==
-        fvc::div(phiHbyA)
-    );
-
-    pEqn.setReference(pRefCell, pRefValue);
-
-    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-    if (pimple.finalNonOrthogonalIter())
-    {
-        phi = phiHbyA - pEqn.flux();
-    }
-}
-
-#include "continuityErrs.H"
-
-p.relax();
-
-U = HbyA - rAU*fvc::grad(p);
-U.correctBoundaryConditions();
-fvOptions.correct(U);
-
-// Correct rhoUf if the mesh is moving
-fvc::correctUf(Uf, U, phi);
-
-// Make the fluxes relative to the mesh motion
-fvc::makeRelative(phi, U);
-

applications/solvers/lagrangian/parcelFoam/parcelFoam.C

@@ -1,154 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    parcelFoam
-
-Group
-    grpLagrangianSolvers
-
-Description
-    Transient solver for incompressible, turbulent flow with kinematic,
-    particle cloud, and surface film modelling.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "dynamicFvMesh.H"
-#include "singlePhaseTransportModel.H"
-#include "turbulentTransportModel.H"
-#include "surfaceFilmModel.H"
-#include "parcelCloudModelList.H"
-#include "fvOptions.H"
-#include "pimpleControl.H"
-#include "CorrectPhi.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::addNote
-    (
-        "Transient solver for incompressible, turbulent flow"
-        " with kinematic particle clouds"
-        " and surface film modelling."
-    );
-
-    #define CREATE_MESH createMeshesPostProcess.H
-    #include "postProcess.H"
-
-    #include "addCheckCaseOptions.H"
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
-    #include "initContinuityErrs.H"
-    #include "createDyMControls.H"
-    #include "createFields.H"
-    #include "createFieldRefs.H"
-    #include "createUfIfPresent.H"
-
-    turbulence->validate();
-
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        #include "readDyMControls.H"
-        #include "CourantNo.H"
-        #include "setMultiRegionDeltaT.H"
-
-        ++runTime;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // Store the particle positions
-        parcels.storeGlobalPositions();
-
-        // Do any mesh changes
-        mesh.update();
-
-        if (pimple.solveFlow() && mesh.changing())
-        {
-            MRF.update();
-
-            if (correctPhi)
-            {
-                // Calculate absolute flux
-                // from the mapped surface velocity
-                phi = mesh.Sf() & Uf();
-
-                #include "../../incompressible/pimpleFoam/correctPhi.H"
-
-                // Make the fluxes relative to the mesh-motion
-                fvc::makeRelative(phi, U);
-            }
-
-            if (checkMeshCourantNo)
-            {
-                #include "meshCourantNo.H"
-            }
-        }
-
-        parcels.evolve();
-        surfaceFilm.evolve();
-
-        if (pimple.solveFlow())
-        {
-            // --- PIMPLE loop
-            while (pimple.loop())
-            {
-                #include "UEqn.H"
-
-                // --- Pressure corrector loop
-                while (pimple.correct())
-                {
-                    #include "pEqn.H"
-                }
-
-                if (pimple.turbCorr())
-                {
-                    laminarTransport.correct();
-                    turbulence->correct();
-                }
-            }
-        }
-
-        runTime.write();
-
-        runTime.printExecutionTime(Info);
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFoam/createClouds.H

@@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-reactingTypeCloud parcels
+basicReactingTypeCloud parcels
 (
     "reactingCloud1",
-    g,
     rho,
     U,
+    g,
     slgThermo
 );

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 reactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
+#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,6 +81,7 @@ int main(int argc, char *argv[])
     #include "createDyMControls.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
+    #include "createRegionControls.H"
     #include "initContinuityErrs.H"
     #include "createRhoUfIfPresent.H"
 
@@ -104,7 +105,7 @@ int main(int argc, char *argv[])
         // so that it can be mapped and used in correctPhi
         // to ensure the corrected phi has the same divergence
         autoPtr<volScalarField> divrhoU;
-        if (pimple.solveFlow() && correctPhi)
+        if (solvePrimaryRegion && correctPhi)
         {
             divrhoU.reset
             (
@@ -132,7 +133,7 @@ int main(int argc, char *argv[])
 
         // Store momentum to set rhoUf for introduced faces.
         autoPtr<volVectorField> rhoU;
-        if (pimple.solveFlow() && rhoUf.valid())
+        if (solvePrimaryRegion && rhoUf.valid())
         {
             rhoU.reset(new volVectorField("rhoU", rho*U));
         }
@@ -143,7 +144,7 @@ int main(int argc, char *argv[])
         // Do any mesh changes
         mesh.update();
 
-        if (pimple.solveFlow() && mesh.changing())
+        if (solvePrimaryRegion && mesh.changing())
         {
             gh = (g & mesh.C()) - ghRef;
             ghf = (g & mesh.Cf()) - ghRef;
@@ -171,7 +172,7 @@ int main(int argc, char *argv[])
         parcels.evolve();
         surfaceFilm.evolve();
 
-        if (pimple.solveFlow())
+        if (solvePrimaryRegion)
         {
             if (pimple.nCorrPIMPLE() <= 1)
             {

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options

@@ -7,6 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H

@@ -1,9 +1,9 @@
 Info<< "\nConstructing " << CLOUD_BASE_TYPE_NAME << " cloud" << endl;
-reactingTypeCloud parcels
+basicReactingTypeCloud parcels
 (
     word(CLOUD_BASE_TYPE_NAME) & "Cloud1",
-    g,
     rho,
     U,
+    g,
     slgThermo
 );

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C

@@ -45,13 +45,13 @@ Description
 #include "cloudMacros.H"
 
 #ifndef CLOUD_BASE_TYPE
-    #define CLOUD_BASE_TYPE reactingMultiphase
+    #define CLOUD_BASE_TYPE ReactingMultiphase
     //#define CLOUD_BASE_TYPE_NAME "reactingMultiphase" Backwards compat
     #define CLOUD_BASE_TYPE_NAME "reacting"
 #endif
 
 #include CLOUD_INCLUDE_FILE(CLOUD_BASE_TYPE)
-#define reactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
+#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/lagrangian/rhoParcelFoam/EEqn.H

@@ -1,41 +0,0 @@
-{
-    volScalarField& he = thermo.he();
-
-    fvScalarMatrix EEqn
-    (
-        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
-      + fvc::ddt(rho, K) + fvc::div(phi, K)
-      + (
-            he.name() == "e"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )
-          : -dpdt
-        )
-      - fvm::laplacian(turbulence->alphaEff(), he)
-     ==
-        rho*(U&g)
-      + parcels.Sh(he)
-      + surfaceFilm.Sh()
-      + radiation->Sh(thermo, he)
-      + Qdot
-      + fvOptions(rho, he)
-    );
-
-    EEqn.relax();
-
-    fvOptions.constrain(EEqn);
-
-    EEqn.solve();
-
-    fvOptions.correct(he);
-
-    thermo.correct();
-    radiation->correct();
-
-    Info<< "T gas min/max   = " << min(T).value() << ", "
-        << max(T).value() << endl;
-}

applications/solvers/lagrangian/rhoParcelFoam/Make/files

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

applications/solvers/lagrangian/rhoParcelFoam/UEqn.H

@@ -1,34 +0,0 @@
-    MRF.correctBoundaryVelocity(U);
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + MRF.DDt(rho, U)
-      + turbulence->divDevRhoReff(U)
-     ==
-        parcels.SU(U)
-      + fvOptions(rho, U)
-    );
-
-    UEqn.relax();
-
-    fvOptions.constrain(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-          ==
-            fvc::reconstruct
-            (
-                (
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                )*mesh.magSf()
-            )
-        );
-
-        fvOptions.correct(U);
-        K = 0.5*magSqr(U);
-    }

applications/solvers/lagrangian/rhoParcelFoam/YEqn.H

@@ -1,51 +0,0 @@
-tmp<fv::convectionScheme<scalar>> mvConvection
-(
-    fv::convectionScheme<scalar>::New
-    (
-        mesh,
-        fields,
-        phi,
-        mesh.divScheme("div(phi,Yi_h)")
-    )
-);
-
-
-{
-    combustion->correct();
-    Qdot = combustion->Qdot();
-    volScalarField Yt(0.0*Y[0]);
-
-    forAll(Y, i)
-    {
-        if (i != inertIndex && composition.active(i))
-        {
-            volScalarField& Yi = Y[i];
-
-            fvScalarMatrix YEqn
-            (
-                fvm::ddt(rho, Yi)
-              + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->muEff(), Yi)
-              ==
-                parcels.SYi(i, Yi)
-              + fvOptions(rho, Yi)
-              + combustion->R(Yi)
-              + surfaceFilm.Srho(i)
-            );
-
-            YEqn.relax();
-
-            fvOptions.constrain(YEqn);
-
-            YEqn.solve(mesh.solver("Yi"));
-
-            fvOptions.correct(Yi);
-
-            Yi.max(0.0);
-            Yt += Yi;
-        }
-    }
-
-    Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
-}

applications/solvers/lagrangian/rhoParcelFoam/createClouds.H

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

applications/solvers/lagrangian/rhoParcelFoam/createFieldRefs.H

@@ -1,5 +0,0 @@
-const volScalarField& T = thermo.T();
-const volScalarField& psi = thermo.psi();
-const label inertIndex(composition.species().find(inertSpecie));
-
-regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();

applications/solvers/lagrangian/rhoParcelFoam/createFields.H

@@ -1,132 +0,0 @@
-#include "createRDeltaT.H"
-
-Info<< "Reading thermophysical properties\n" << endl;
-autoPtr<rhoReactionThermo> pThermo(rhoReactionThermo::New(mesh));
-rhoReactionThermo& thermo = pThermo();
-thermo.validate(args.executable(), "h", "e");
-
-SLGThermo slgThermo(mesh, thermo);
-
-basicSpecieMixture& composition = thermo.composition();
-PtrList<volScalarField>& Y = composition.Y();
-
-const word inertSpecie(thermo.get<word>("inertSpecie"));
-if
-(
-    !composition.species().found(inertSpecie)
-  && composition.species().size() > 0
-)
-{
-    FatalIOErrorIn(args.executable().c_str(), thermo)
-        << "Inert specie " << inertSpecie << " not found in available species "
-        << composition.species()
-        << exit(FatalIOError);
-}
-
-Info<< "Creating field rho\n" << endl;
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
-
-volScalarField& p = thermo.p();
-
-Info<< "\nReading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "compressibleCreatePhi.H"
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-Info<< "Creating combustion model\n" << endl;
-autoPtr<CombustionModel<rhoReactionThermo>> combustion
-(
-    CombustionModel<rhoReactionThermo>::New(thermo, turbulence())
-);
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-// Force p_rgh to be consistent with p
-p_rgh = p - rho*gh;
-
-pressureControl pressureControl(p, rho, pimple.dict(), false);
-
-mesh.setFluxRequired(p_rgh.name());
-
-Info<< "Creating multi-variate interpolation scheme\n" << endl;
-multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-forAll(Y, i)
-{
-    fields.add(Y[i]);
-}
-fields.add(thermo.he());
-
-volScalarField Qdot
-(
-    IOobject
-    (
-        "Qdot",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar(dimEnergy/dimVolume/dimTime, Zero)
-);
-
-#include "createDpdt.H"
-
-#include "createK.H"
-
-#include "createMRF.H"
-#include "createRadiationModel.H"
-#include "createClouds.H"
-#include "createSurfaceFilmModel.H"
-#include "createFvOptions.H"

applications/solvers/lagrangian/rhoParcelFoam/createMeshesPostProcess.H

@@ -1,35 +0,0 @@
-#include "createMesh.H"
-
-dictionary filmDict;
-
-IOobject io
-(
-    "surfaceFilmProperties",
-    mesh.time().constant(),
-    mesh,
-    IOobject::READ_IF_PRESENT,
-    IOobject::NO_WRITE,
-    false
-);
-
-if (io.typeHeaderOk<IOdictionary>())
-{
-    IOdictionary propDict(io);
-
-    filmDict = std::move(propDict);
-
-    const word filmRegionName = filmDict.get<word>("region");
-
-    fvMesh filmMesh
-    (
-        IOobject
-        (
-            filmRegionName,
-            runTime.timeName(),
-            runTime,
-            IOobject::MUST_READ
-        )
-    );
-}
-
-

applications/solvers/lagrangian/rhoParcelFoam/createSurfaceFilmModel.H

@@ -1,6 +0,0 @@
-Info<< "\nConstructing surface film model" << endl;
-
-autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
-(
-    regionModels::surfaceFilmModel::New(mesh, g)
-);

applications/solvers/lagrangian/rhoParcelFoam/pEqn.H

@@ -1,96 +0,0 @@
-if (!pimple.SIMPLErho())
-{
-    rho = thermo.rho();
-}
-
-// Thermodynamic density needs to be updated by psi*d(p) after the
-// pressure solution
-const volScalarField psip0(psi*p);
-
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    (
-        fvc::flux(rho*HbyA)
-      + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
-    )
-  + phig
-);
-
-fvc::makeRelative(phiHbyA, rho, U);
-MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-// Update the pressure BCs to ensure flux consistency
-constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-fvScalarMatrix p_rghDDtEqn
-(
-    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-  + fvc::div(phiHbyA)
- ==
-    parcels.Srho()
-  + surfaceFilm.Srho()
-  + fvOptions(psi, p_rgh, rho.name())
-);
-
-while (pimple.correctNonOrthogonal())
-{
-    fvScalarMatrix p_rghEqn
-    (
-        p_rghDDtEqn
-      - fvm::laplacian(rhorAUf, p_rgh)
-    );
-
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-    if (pimple.finalNonOrthogonalIter())
-    {
-        phi = phiHbyA + p_rghEqn.flux();
-
-        // Explicitly relax pressure for momentum corrector
-        p_rgh.relax();
-
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
-        U.correctBoundaryConditions();
-        fvOptions.correct(U);
-        K = 0.5*magSqr(U);
-    }
-}
-
-p = p_rgh + rho*gh;
-
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-if (pressureControl.limit(p))
-{
-    p.correctBoundaryConditions();
-    rho = thermo.rho();
-    p_rgh = p - rho*gh;
-}
-else if (pimple.SIMPLErho())
-{
-    rho = thermo.rho();
-}
-
-// Correct rhoUf if the mesh is moving
-fvc::correctRhoUf(rhoUf, rho, U, phi);
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-
-    if (mesh.moving())
-    {
-        dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
-    }
-}

applications/solvers/lagrangian/rhoParcelFoam/rhoEqn.H

@@ -1,50 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2017 OpenFOAM Foundation
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Global
-    rhoEqn
-
-Description
-    Solve the continuity for density.
-
-\*---------------------------------------------------------------------------*/
-
-{
-    fvScalarMatrix rhoEqn
-    (
-        fvm::ddt(rho)
-      + fvc::div(phi)
-      ==
-        parcels.Srho(rho)
-      + surfaceFilm.Srho()
-      + fvOptions(rho)
-    );
-
-    rhoEqn.solve();
-
-    fvOptions.correct(rho);
-}
-
-// ************************************************************************* //

applications/solvers/lagrangian/rhoParcelFoam/rhoParcelFoam.C

@@ -1,206 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2020 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    rhoParcelFoam
-
-Group
-    grpLagrangianSolvers
-
-Description
-    Transient solver for compressible, turbulent flow with a reacting,
-    multiphase particle cloud, and surface film modelling.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "dynamicFvMesh.H"
-#include "turbulentFluidThermoModel.H"
-#include "surfaceFilmModel.H"
-#include "rhoReactionThermo.H"
-#include "CombustionModel.H"
-#include "radiationModel.H"
-#include "SLGThermo.H"
-#include "fvOptions.H"
-#include "pimpleControl.H"
-#include "pressureControl.H"
-#include "CorrectPhi.H"
-#include "localEulerDdtScheme.H"
-#include "fvcSmooth.H"
-
-#include "parcelCloudModelList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::addNote
-    (
-        "Transient solver for compressible, turbulent flow"
-        " with lagrangian parcel clouds and surface film modelling."
-    );
-
-    #define CREATE_MESH createMeshesPostProcess.H
-    #include "postProcess.H"
-
-    #include "addCheckCaseOptions.H"
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
-    #include "createDyMControls.H"
-    #include "createFields.H"
-    #include "createFieldRefs.H"
-    #include "initContinuityErrs.H"
-    #include "createRhoUfIfPresent.H"
-
-    turbulence->validate();
-
-    if (!LTS)
-    {
-        #include "compressibleCourantNo.H"
-        #include "setInitialDeltaT.H"
-    }
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        #include "readDyMControls.H"
-
-        // Store divrhoU from the previous mesh
-        // so that it can be mapped and used in correctPhi
-        // to ensure the corrected phi has the same divergence
-        autoPtr<volScalarField> divrhoU;
-        if (pimple.solveFlow() && correctPhi)
-        {
-            divrhoU.reset
-            (
-                new volScalarField
-                (
-                    "divrhoU",
-                    fvc::div(fvc::absolute(phi, rho, U))
-                )
-            );
-        }
-
-        if (LTS)
-        {
-            #include "setRDeltaT.H"
-        }
-        else
-        {
-            #include "compressibleCourantNo.H"
-            #include "setMultiRegionDeltaT.H"
-        }
-
-        ++runTime;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // Store momentum to set rhoUf for introduced faces.
-        autoPtr<volVectorField> rhoU;
-        if (pimple.solveFlow() && rhoUf.valid())
-        {
-            rhoU.reset(new volVectorField("rhoU", rho*U));
-        }
-
-        // Store the particle positions
-        parcels.storeGlobalPositions();
-
-        // Do any mesh changes
-        mesh.update();
-
-        if (pimple.solveFlow() && mesh.changing())
-        {
-            gh = (g & mesh.C()) - ghRef;
-            ghf = (g & mesh.Cf()) - ghRef;
-
-            MRF.update();
-
-            if (correctPhi)
-            {
-                // Calculate absolute flux
-                // from the mapped surface velocity
-                phi = mesh.Sf() & rhoUf();
-
-                #include "../../compressible/rhoPimpleFoam/correctPhi.H"
-
-                // Make the fluxes relative to the mesh-motion
-                fvc::makeRelative(phi, rho, U);
-            }
-
-            if (checkMeshCourantNo)
-            {
-                #include "meshCourantNo.H"
-            }
-        }
-
-        parcels.evolve();
-        surfaceFilm.evolve();
-
-        if (pimple.solveFlow())
-        {
-            if (pimple.nCorrPIMPLE() <= 1)
-            {
-                #include "rhoEqn.H"
-            }
-
-            // --- PIMPLE loop
-            while (pimple.loop())
-            {
-                #include "UEqn.H"
-                #include "YEqn.H"
-                #include "EEqn.H"
-
-                // --- Pressure corrector loop
-                while (pimple.correct())
-                {
-                    #include "pEqn.H"
-                }
-
-                if (pimple.turbCorr())
-                {
-                    turbulence->correct();
-                }
-            }
-
-            rho = thermo.rho();
-        }
-
-        runTime.write();
-
-        runTime.printExecutionTime(Info);
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/lagrangian/rhoParcelFoam/setMultiRegionDeltaT.H

@@ -1,55 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2017 OpenFOAM Foundation
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Global
-    setMultiRegionDeltaT
-
-Description
-    Reset the timestep to maintain a constant maximum Courant numbers.
-    Reduction of time-step is immediate, but increase is damped to avoid
-    unstable oscillations.
-
-\*---------------------------------------------------------------------------*/
-
-if (adjustTimeStep)
-{
-    const scalar maxDeltaTFact =
-        min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
-    const scalar deltaTFact =
-        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
-
-    runTime.setDeltaT
-    (
-        min
-        (
-            deltaTFact*runTime.deltaTValue(),
-            maxDeltaT
-        )
-    );
-
-    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
-}
-
-// ************************************************************************* //

applications/solvers/lagrangian/rhoParcelFoam/setRDeltaT.H

@@ -1,137 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-{
-    volScalarField& rDeltaT = trDeltaT.ref();
-
-    const dictionary& pimpleDict = pimple.dict();
-
-    // Maximum flow Courant number
-    scalar maxCo(pimpleDict.get<scalar>("maxCo"));
-
-    // Maximum time scale
-    scalar maxDeltaT(pimpleDict.getOrDefault<scalar>("maxDeltaT", GREAT));
-
-    // Smoothing parameter (0-1) when smoothing iterations > 0
-    scalar rDeltaTSmoothingCoeff
-    (
-        pimpleDict.getOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
-    );
-
-    // Damping coefficient (1-0)
-    scalar rDeltaTDampingCoeff
-    (
-        pimpleDict.getOrDefault<scalar>("rDeltaTDampingCoeff", 0.2)
-    );
-
-    // Maximum change in cell temperature per iteration
-    // (relative to previous value)
-    scalar alphaTemp(pimpleDict.getOrDefault("alphaTemp", 0.05));
-
-
-    Info<< "Time scales min/max:" << endl;
-
-    // Cache old reciprocal time scale field
-    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
-
-    // Flow time scale
-    {
-        rDeltaT.ref() =
-        (
-            fvc::surfaceSum(mag(phi))()()
-           /((2*maxCo)*mesh.V()*rho())
-        );
-
-        // Limit the largest time scale
-        rDeltaT.max(1/maxDeltaT);
-
-        Info<< "    Flow        = "
-            << gMin(1/rDeltaT.primitiveField()) << ", "
-            << gMax(1/rDeltaT.primitiveField()) << endl;
-    }
-
-    // Reaction source time scale
-    {
-        volScalarField::Internal rDeltaTT
-        (
-            mag
-            (
-                parcels.hsTrans()/(mesh.V()*runTime.deltaT())
-              + Qdot()
-            )
-           /(
-               alphaTemp
-              *rho()
-              *thermo.Cp()()()
-              *T()
-           )
-        );
-
-        Info<< "    Temperature = "
-            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
-            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
-
-        rDeltaT.ref() = max
-        (
-            rDeltaT(),
-            rDeltaTT
-        );
-    }
-
-    // Update the boundary values of the reciprocal time-step
-    rDeltaT.correctBoundaryConditions();
-
-    // Spatially smooth the time scale field
-    if (rDeltaTSmoothingCoeff < 1.0)
-    {
-        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
-    }
-
-    // Limit rate of change of time scale
-    // - reduce as much as required
-    // - only increase at a fraction of old time scale
-    if
-    (
-        rDeltaTDampingCoeff < 1.0
-     && runTime.timeIndex() > runTime.startTimeIndex() + 1
-    )
-    {
-        rDeltaT = max
-        (
-            rDeltaT,
-            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
-        );
-    }
-
-    Info<< "    Overall     = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
-}
-
-
-// ************************************************************************* //

applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options

@@ -7,6 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@@ -46,6 +47,7 @@ EXE_LIBS = \
     -lsurfaceFilmModels \
     -lcombustionModels \
     -lsampling \
+    -lcoalCombustion \
     -lregionFaModels \
     -lfiniteArea \
     -lfaOptions

applications/solvers/lagrangian/simpleCoalParcelFoam/createClouds.H

@@ -1,9 +1,9 @@
 Info<< "\nConstructing coal cloud" << endl;
-reactingMultiphaseCloud parcels
+coalCloud parcels
 (
     "reactingCloud1",
-    g,
     rho,
     U,
+    g,
     slgThermo
 );

applications/solvers/lagrangian/simpleCoalParcelFoam/simpleCoalParcelFoam.C

@@ -37,7 +37,7 @@ Description
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "reactingMultiphaseCloud.H"
+#include "coalCloud.H"
 #include "rhoReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"

applications/solvers/lagrangian/sprayFoam/Make/options

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

applications/solvers/lagrangian/sprayFoam/createClouds.H

@@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-sprayCloud parcels
+basicSprayCloud parcels
 (
     "sprayCloud",
-    g,
     rho,
     U,
+    g,
     slgThermo
 );

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

@@ -39,6 +39,7 @@ 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 "sprayCloud.H"
+#include "basicSprayCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"

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

@@ -34,6 +34,7 @@ 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,6 +40,7 @@ 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 "sprayCloud.H"
+#include "basicSprayCloud.H"
 #include "psiReactionThermo.H"
 #include "CombustionModel.H"
 #include "radiationModel.H"