TUT: added two unsteady adjoint optimisation tutorials

showcasing the use of compressedFullStorage and binomialCheckPointing

.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=230110
+api=2206
+patch=0

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/acoustic/acousticFoam/Make/options

@@ -11,5 +11,4 @@ EXE_LIBS = \
     -lfvOptions \
     -lmeshTools \
     -lsampling \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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/PDRFoam/PDRFoamAutoRefine.C

@@ -116,7 +116,7 @@ int main(int argc, char *argv[])
             "normalisedGradP",
             tmagGradP()/max(tmagGradP())
         );
-        normalisedGradP.writeOpt(IOobject::AUTO_WRITE);
+        normalisedGradP.writeOpt() = IOobject::AUTO_WRITE;
         tmagGradP.clear();
 
         ++runTime;

applications/solvers/combustion/fireFoam/Make/options

@@ -1,8 +1,8 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \

applications/solvers/compressible/rhoPimpleFoam/Make/options

@@ -9,7 +9,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -26,5 +26,4 @@ EXE_LIBS = \
     -ltopoChangerFvMesh \
     -lsampling \
     -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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/compressible/rhoSimpleFoam/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -21,5 +21,4 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModels \
     -lthermoTools \
     -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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/Make/options

@@ -10,7 +10,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -28,5 +28,4 @@ EXE_LIBS = \
     -ldynamicMesh \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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/buoyantSimpleFoam/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -22,5 +22,4 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModels \
     -lthermoTools \
     -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -19,7 +19,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 
 EXE_LIBS = \
@@ -40,5 +40,4 @@ EXE_LIBS = \
     -lthermoTools \
     -lradiationModels \
     -lregionModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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/Make/options

@@ -36,5 +36,5 @@ EXE_LIBS = \
     -lfaOptions \
     -lregionModels \
     -lsampling \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
+

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/Make/options

@@ -8,7 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/regionFaModels\lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -22,5 +22,4 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
     -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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/coalChemistryFoam/Make/options

@@ -1,7 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \

applications/solvers/lagrangian/reactingParcelFoam/Make/options

@@ -2,10 +2,10 @@ EXE_INC = \
     -I../reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \

applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options

@@ -2,8 +2,8 @@ EXE_INC = \
     -I.. \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

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

@@ -1,7 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \

applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options

@@ -2,7 +2,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \

applications/solvers/lagrangian/sprayFoam/Make/options

@@ -2,8 +2,8 @@ EXE_INC = \
     -I../reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \

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

@@ -4,8 +4,8 @@ EXE_INC = \
     -I../../reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \

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

@@ -2,7 +2,7 @@ EXE_INC = \
     -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam/simpleReactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \

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

@@ -4,8 +4,8 @@ EXE_INC = \
     -I../../../compressible/rhoPimpleFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \

applications/solvers/multiphase/MPPICInterFoam/Make/options

@@ -42,5 +42,4 @@ EXE_LIBS = \
     -lregionModels \
     -lregionFaModels \
     -lsurfaceFilmModels \
-    -lfiniteArea \
     -lcompressibleTwoPhaseMixtureTurbulenceModels

applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H

@@ -129,5 +129,10 @@ compressibleInterPhaseTransportModel turbulence
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
 
+surfaceScalarField faceMask
+(
+    localMin<scalar>(mesh).interpolate(cellMask)
+);
+
 // Create bool field with interpolated cells
 #include "createInterpolatedCells.H"

applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/overCompressibleInterDyMFoam.C

@@ -80,6 +80,9 @@ int main(int argc, char *argv[])
     #include "setRootCaseLists.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
+    pimpleControl pimple(mesh);
+
+    #include "createTimeControls.H"
     #include "createDyMControls.H"
     #include "createFields.H"
 
@@ -106,7 +109,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readDyMControls.H"
+        #include "readControls.H"
 
         if (LTS)
         {
@@ -151,10 +154,40 @@ int main(int argc, char *argv[])
                     // Update cellMask field for blocking out hole cells
                     #include "setCellMask.H"
                     #include "setInterpolatedCells.H"
-                    #include "correctPhiFaceMask.H"
+
+                    faceMask =
+                        localMin<scalar>(mesh).interpolate(cellMask.oldTime());
+
+                    // Zero Uf on old faceMask (H-I)
+                    Uf *= faceMask;
+
+                    const surfaceVectorField Uint(fvc::interpolate(U));
+                    // Update Uf and phi on new C-I faces
+                    Uf += (1-faceMask)*Uint;
+
+                    // Update Uf boundary
+                    forAll(Uf.boundaryField(), patchI)
+                    {
+                        Uf.boundaryFieldRef()[patchI] =
+                            Uint.boundaryField()[patchI];
+                    }
+
+                    phi = mesh.Sf() & Uf;
+
+                    // Correct phi on individual regions
+                    if (correctPhi)
+                    {
+                         #include "correctPhi.H"
+                    }
 
                     mixture.correct();
 
+                    // Zero phi on current H-I
+                    faceMask = localMin<scalar>(mesh).interpolate(cellMask);
+
+                    phi *= faceMask;
+                    U *= cellMask;
+
                     // Make the flux relative to the mesh motion
                     fvc::makeRelative(phi, U);
 
@@ -169,6 +202,10 @@ int main(int argc, char *argv[])
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
 
+            const surfaceScalarField faceMask
+            (
+                localMin<scalar>(mesh).interpolate(cellMask)
+            );
             rhoPhi *= faceMask;
 
             turbulence.correctPhasePhi();

applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H

@@ -10,6 +10,19 @@
         fvc::flux(HbyA)
     );
 
+    if (ddtCorr)
+    {
+        surfaceScalarField faceMaskOld
+        (
+            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+        );
+        phiHbyA +=
+            MRF.zeroFilter
+            (
+                fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf)
+            );
+    }
+
     MRF.makeRelative(phiHbyA);
 
     surfaceScalarField phig

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleIO.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,7 +59,7 @@ Foam::DTRMParticle::DTRMParticle
 {
     if (readFields)
     {
-        if (is.format() == IOstreamOption::ASCII)
+        if (is.format() == IOstream::ASCII)
         {
             is >> p0_ >> p1_ >> I0_ >> I_ >> dA_ >> transmissiveId_;
         }
@@ -115,7 +115,7 @@ void Foam::DTRMParticle::writeProperties
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const DTRMParticle& p)
 {
-    if (os.format() == IOstreamOption::ASCII)
+    if (os.format() == IOstream::ASCII)
     {
         os  << static_cast<const particle&>(p)
             << token::SPACE << p.p0_

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C

@@ -157,7 +157,12 @@ void Foam::radiation::laserDTRM::initialiseReflection()
             );
         }
 
-        reflectionSwitch_ = returnReduceOr(reflections_.size());
+        if (reflections_.size())
+        {
+            reflectionSwitch_ = true;
+        }
+
+        reflectionSwitch_ = returnReduce(reflectionSwitch_, orOp<bool>());
     }
 }
 
@@ -294,12 +299,14 @@ void Foam::radiation::laserDTRM::initialise()
                     DTRMCloud_.addParticle(pPtr);
                 }
 
-                if (nMissed < 10 && returnReduceAnd(cellI < 0))
+                if (returnReduce(cellI, maxOp<label>()) == -1)
                 {
-                    ++nMissed;
-                    WarningInFunction
-                        << "Cannot find owner cell for focalPoint at "
-                        << p0 << endl;
+                    if (++nMissed <= 10)
+                    {
+                        WarningInFunction
+                            << "Cannot find owner cell for focalPoint at "
+                            << p0 << endl;
+                    }
                 }
             }
         }
@@ -692,7 +699,7 @@ void Foam::radiation::laserDTRM::calculate()
 
             for (label pointi = 0; pointi < lines.size(); pointi += 2)
             {
-                os.writeLine(lines[pointi], lines[pointi+1]);
+                os.write(linePointRef(lines[pointi], lines[pointi+1]));
             }
         }
 

applications/solvers/multiphase/interFoam/overInterDyMFoam/CourantNo.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2011 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,10 +36,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
     (

applications/solvers/multiphase/interFoam/overInterDyMFoam/alphaSuSp.H

@@ -7,5 +7,3 @@ volScalarField::Internal divU
   ? fvc::div(phiCN() + mesh.phi())
   : fvc::div(phiCN())
 );
-
-divU *= interpolatedCells*cellMask;

applications/solvers/multiphase/interFoam/overInterDyMFoam/createDyMControls.H

@@ -0,0 +1,30 @@
+bool correctPhi
+(
+    pimple.dict().getOrDefault("correctPhi", true)
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().getOrDefault("checkMeshCourantNo", false)
+);
+
+bool moveMeshOuterCorrectors
+(
+    pimple.dict().getOrDefault("moveMeshOuterCorrectors", false)
+);
+
+
+bool massFluxInterpolation
+(
+    pimple.dict().getOrDefault("massFluxInterpolation", false)
+);
+
+bool adjustFringe
+(
+    pimple.dict().getOrDefault("oversetAdjustPhi", false)
+);
+
+bool ddtCorr
+(
+    pimple.dict().getOrDefault("ddtCorr", true)
+);

applications/solvers/multiphase/interFoam/overInterDyMFoam/createFields.H

@@ -40,12 +40,8 @@ volVectorField U
     nonInt.insert("HbyA");
     nonInt.insert("grad(p_rgh)");
     nonInt.insert("nHat");
-    nonInt.insert("surfaceIntegrate(nHatf)");
-    nonInt.insert("surfaceIntegrate(phi+meshPhi)");
     nonInt.insert("surfaceIntegrate(phi)");
     nonInt.insert("surfaceIntegrate(phiHbyA)");
-    nonInt.insert("surfaceSum(((S|magSf)*S)");
-    nonInt.insert("surfaceIntegrate(((rAUf*magSf)*snGradCorr(p_rgh)))");
     nonInt.insert("cellMask");
     nonInt.insert("cellDisplacement");
     nonInt.insert("interpolatedCells");

applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2014 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,11 +49,14 @@ Description
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
+#include "CorrectPhi.H"
 #include "fvcSmooth.H"
 #include "cellCellStencilObject.H"
 #include "localMin.H"
+#include "interpolationCellPoint.H"
+#include "transform.H"
+#include "fvMeshSubset.H"
 #include "oversetAdjustPhi.H"
-#include "oversetPatchPhiErr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +76,8 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "initContinuityErrs.H"
-
+    pimpleControl pimple(mesh);
+    #include "createTimeControls.H"
     #include "createDyMControls.H"
     #include "createFields.H"
     #include "createAlphaFluxes.H"
@@ -93,8 +97,11 @@ int main(int argc, char *argv[])
         dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
     );
 
+    if (correctPhi)
+    {
+        #include "correctPhi.H"
+    }
     #include "createUf.H"
-    #include "createControls.H"
 
     #include "setCellMask.H"
     #include "setInterpolatedCells.H"
@@ -112,8 +119,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readDyMControls.H"
-        #include "readOversetDyMControls.H"
+        #include "readControls.H"
 
         if (LTS)
         {
@@ -152,17 +158,50 @@ int main(int argc, char *argv[])
                         talphaPhi1Corr0.clear();
                     }
 
-                    // Update cellMask field for blocking out hole cells
-                    #include "setCellMask.H"
-                    #include "setInterpolatedCells.H"
-                    #include "correctPhiFaceMask.H"
-
                     gh = (g & mesh.C()) - ghRef;
                     ghf = (g & mesh.Cf()) - ghRef;
 
+                    // Update cellMask field for blocking out hole cells
+                    #include "setCellMask.H"
+                    #include "setInterpolatedCells.H"
+
+                    const surfaceScalarField faceMaskOld
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+                    );
+
+                    // Zero Uf on old faceMask (H-I)
+                    Uf *= faceMaskOld;
+
+                    const surfaceVectorField Uint(fvc::interpolate(U));
+                    // Update Uf and phi on new C-I faces
+                    Uf += (1-faceMaskOld)*Uint;
+
+                    // Update Uf boundary
+                    forAll(Uf.boundaryField(), patchI)
+                    {
+                        Uf.boundaryFieldRef()[patchI] =
+                            Uint.boundaryField()[patchI];
+                    }
+
+                    phi = mesh.Sf() & Uf;
+
+                    // Correct phi on individual regions
+                    if (correctPhi)
+                    {
+                         #include "correctPhi.H"
+                    }
 
                     mixture.correct();
 
+                    // Zero phi on current H-I
+                    const surfaceScalarField faceMask
+                    (
+                        localMin<scalar>(mesh).interpolate(cellMask)
+                    );
+                    phi *= faceMask;
+                    U   *= cellMask;
+
                     // Make the flux relative to the mesh motion
                     fvc::makeRelative(phi, U);
 
@@ -174,14 +213,14 @@ int main(int argc, char *argv[])
                 }
             }
 
-            if (adjustFringe)
-            {
-                oversetAdjustPhi(phi, U, zoneIdMass);
-            }
 
             #include "alphaControls.H"
             #include "alphaEqnSubCycle.H"
 
+            const surfaceScalarField faceMask
+            (
+                localMin<scalar>(mesh).interpolate(cellMask)
+            );
             rhoPhi *= faceMask;
 
             mixture.correct();

applications/solvers/multiphase/interFoam/overInterDyMFoam/pEqn.H

@@ -1,31 +1,64 @@
 {
     rAU = 1.0/UEqn.A();
+    //mesh.interpolate(rAU);
+
+    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField H("H", UEqn.H());
 
     volVectorField HbyA("HbyA", U);
-
+    //HbyA = rAU*UEqn.H();
     HbyA = constrainHbyA(rAU*H, U, p_rgh);
 
+    if (massFluxInterpolation)
+    {
+        #include "interpolatedFaces.H"
+    }
+
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 
+    if (ddtCorr)
+    {
+        surfaceScalarField faceMaskOld
+        (
+            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+        );
+        phiHbyA +=
+            fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf);
+    }
     MRF.makeRelative(phiHbyA);
 
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
+
+    if (adjustFringe)
+    {
+        fvc::makeRelative(phiHbyA, U);
+        oversetAdjustPhi(phiHbyA, U);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
+
     surfaceScalarField phig
     (
         (
             mixture.surfaceTensionForce()
-          - ghf*fvc::snGrad(cellMask*rho)
-        )*rAUf*faceMask*mesh.magSf()
+          - ghf*fvc::snGrad(rho)
+        )*faceMask*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
     if (adjustFringe)
     {
-        oversetAdjustPhi(phiHbyA, U, zoneIdMass);
+        fvc::makeRelative(phiHbyA, U);
+        oversetAdjustPhi(phiHbyA, U);
+        fvc::makeAbsolute(phiHbyA, U);
     }
 
     // Update the pressure BCs to ensure flux consistency
@@ -57,10 +90,6 @@
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
-        if (oversetPatchErrOutput)
-        {
-            oversetPatchPhiErr(p_rghEqn, phiHbyA);
-        }
     }
 
     #include "continuityErrs.H"

applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H

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

applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/correctPhi.H

@@ -0,0 +1,11 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    surfaceScalarField("rAUf", fvc::interpolate(rAU)),
+    divU,
+    pimple
+);
+
+#include "continuityErrs.H"

applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H

@@ -149,5 +149,10 @@ surfaceScalarField alphaPhi10
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
 
+surfaceScalarField faceMask
+(
+    localMin<scalar>(mesh).interpolate(cellMask)
+);
+
 // Create bool field with interpolated cells
 #include "createInterpolatedCells.H"

applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/overInterPhaseChangeDyMFoam.C

@@ -82,7 +82,9 @@ int main(int argc, char *argv[])
     #include "setRootCaseLists.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
+    pimpleControl pimple(mesh);
 
+    #include "createTimeControls.H"
     #include "createDyMControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -116,7 +118,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readDyMControls.H"
+        #include "readControls.H"
 
         // Store divU from the previous mesh so that it can be mapped
         // and used in correctPhi to ensure the corrected phi has the
@@ -151,10 +153,40 @@ int main(int argc, char *argv[])
                      // Update cellMask field for blocking out hole cells
                     #include "setCellMask.H"
                     #include "setInterpolatedCells.H"
-                    #include "correctPhiFaceMask.H"
+
+                    faceMask =
+                        localMin<scalar>(mesh).interpolate(cellMask.oldTime());
+
+
+                    // Zero Uf on old faceMask (H-I)
+                    Uf *= faceMask;
+
+                    const surfaceVectorField Uint(fvc::interpolate(U));
+                    // Update Uf and phi on new C-I faces
+                    Uf += (1-faceMask)*Uint;
+
+                    // Update Uf boundary
+                    forAll(Uf.boundaryField(), patchI)
+                    {
+                        Uf.boundaryFieldRef()[patchI] =
+                            Uint.boundaryField()[patchI];
+                    }
+
+                    phi = mesh.Sf() & Uf;
+
+                    if (correctPhi)
+                    {
+                        #include "correctPhi.H"
+                    }
 
                     mixture->correct();
 
+                    // Zero phi on current H-I
+                    faceMask = localMin<scalar>(mesh).interpolate(cellMask);
+
+                    phi *= faceMask;
+                    U   *= cellMask;
+
                     // Make the flux relative to the mesh motion
                     fvc::makeRelative(phi, U);
                 }
@@ -182,6 +214,10 @@ int main(int argc, char *argv[])
             mixture->correct();
 
             #include "alphaEqnSubCycle.H"
+            const surfaceScalarField faceMask
+            (
+                localMin<scalar>(mesh).interpolate(cellMask)
+            );
             rhoPhi *= faceMask;
 
             interface.correct();

applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H

@@ -8,6 +8,16 @@
         fvc::flux(HbyA)
     );
 
+
+    if (ddtCorr)
+    {
+        surfaceScalarField faceMaskOld
+        (
+            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+        );
+        phiHbyA += faceMaskOld*fvc::ddtCorr(U, Uf);
+    }
+
     if (p_rgh.needReference())
     {
         fvc::makeRelative(phiHbyA, U);

applications/test/00-dummy/Allwmake

@@ -7,7 +7,6 @@ cd "${0%/*}" || exit                                # Run from this directory
 
 library/Allwmake
 
-# Does not use libOpenFOAM or libPstream...
-wmake -no-openfoam
+wmake
 
 #------------------------------------------------------------------------------

applications/test/00-dummy/Make/files

@@ -1,4 +1,4 @@
 /* Install into FOAM_APPBIN to support 'fake' builds (to test packaging) */
-Test-dummyLib.cpp
+Test-dummyLib.C
 
 EXE = $(FOAM_APPBIN)/Test-dummyLib

applications/test/00-dummy/library/Make/files

@@ -1,4 +1,4 @@
 /* Install into central FOAM_LIBBIN for 'fake' builds (test packaging) */
-dummyLib.cpp
+dummyLib.C
 
 LIB = $(FOAM_LIBBIN)/libTestDummy-OpenFOAM

applications/test/00-dummy/library/dummy/Make/files

@@ -1,4 +1,4 @@
 /* Install into central FOAM_LIBBIN/dummy for 'fake' builds (test packaging) */
-dummyMpiLib.cpp
+dummyMpiLib.C
 
 LIB = $(FOAM_LIBBIN)/dummy/libTestDummyMpi

applications/test/00-dummy/library/mpi/Make/files

@@ -1,4 +1,4 @@
 /* Install into central FOAM_MPI_LIBBIN for 'fake' builds (test packaging) */
-dummyMpiLib.cpp
+dummyMpiLib.C
 
 LIB = $(FOAM_MPI_LIBBIN)/libTestDummyMpi

applications/test/00-machine-sizes/Allwmake

@@ -1,8 +0,0 @@
-#!/bin/sh
-cd "${0%/*}" || exit                                # Run from this directory
-#------------------------------------------------------------------------------
-
-# Does not use libOpenFOAM or libPstream...
-wmake -no-openfoam
-
-#------------------------------------------------------------------------------

applications/test/00-machine-sizes/Test-machine-sizes.cpp

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,8 +59,7 @@ void print(const char* name, bool showLimits = true)
 
 int main(int argc, char *argv[])
 {
-    std::cout<< "c++ = " << __cplusplus << '\n';
-    std::cout<< "machine sizes\n---\n\n";
+    std::cout<<"machine sizes\n---\n\n";
 
     print<short>("short");
     print<int>("int");

applications/test/00-openmp/Allwmake

@@ -1,8 +0,0 @@
-#!/bin/sh
-cd "${0%/*}" || exit                                # Run from this directory
-#------------------------------------------------------------------------------
-
-# Does not use libOpenFOAM or libPstream...
-wmake -no-openfoam
-
-#------------------------------------------------------------------------------

applications/test/00-openmp/Make/options

@@ -1,4 +0,0 @@
-EXE_INC = $(COMP_OPENMP)
-
-/* Mostly do not need to explicitly link openmp libraries */
-/* EXE_LIBS = $(LINK_OPENMP) */

applications/test/CircularBuffer/Test-CircularBuffer.C

@@ -64,23 +64,23 @@ inline Ostream& report
 int main(int argc, char *argv[])
 {
     CircularBuffer<label> buf1(1); report(buf1);
-    buf1.push_back(10); report(buf1);
+    buf1.append(10); report(buf1);
 
     Info<< buf1.range_one() << nl;
 
-    buf1.push_back(20); report(buf1);
-    buf1.push_back(30); report(buf1);
+    buf1.append(20); report(buf1);
+    buf1.append(30); report(buf1);
     buf1.push_back(40); report(buf1);
     buf1.push_front(-50); report(buf1);
-    buf1.push_back(60); report(buf1);
-    buf1.push_back(labelList({70,80,90})); report(buf1);
+    buf1.append(60); report(buf1);
+    buf1.append(labelList({70,80,90})); report(buf1);
 
     Info<< nl << "access: " << buf1 << nl;
 
     Info<< buf1[-12] << nl;
 
     Info<< "found: " << buf1.found(40) << nl;
-    buf1.push_uniq(100); report(buf1);
+    buf1.appendUniq(100); report(buf1);
 
     buf1 = Zero; report(buf1);
 
@@ -88,11 +88,11 @@ int main(int argc, char *argv[])
 
     while (buf1.size() > 2)
     {
-        buf1.pop_front();
+        (void) buf1.pop_front();
     }
     report(buf1);
 
-    buf1.push_back(identity(5)); report(buf1);
+    buf1.append(identity(5)); report(buf1);
 
     buf1.info(Info);
     Info<< buf1 << nl;
@@ -123,8 +123,6 @@ int main(int argc, char *argv[])
         Info<< endl;
     }
 
-    Info<< nl << "list: " << flatOutput(buf2.list()) << nl;
-
     Info<< "normal: " << flatOutput(buf2) << nl;
     buf2.reverse();
     Info<< "reverse: " << flatOutput(buf2) << nl;

applications/test/CompactIOList/Test-CompactIOList.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,8 +49,8 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
 
-    IOstreamOption streamOpt(IOstreamOption::BINARY);
-    // IOstreamOption streamOpt(IOstreamOption::ASCII);
+    IOstream::streamFormat format = IOstream::BINARY;
+    // IOstream::streamFormat format = IOstream::ASCII;
 
     const label size = 20000000;
 
@@ -85,7 +85,11 @@ int main(int argc, char *argv[])
             << runTime.cpuTimeIncrement() << " s" << nl << endl;
 
 
-        faces2.writeObject(streamOpt, true);
+        faces2.writeObject
+        (
+            IOstreamOption(format),
+            true
+        );
 
         Info<< "Written old format faceList in = "
             << runTime.cpuTimeIncrement() << " s" << nl << endl;
@@ -141,7 +145,11 @@ int main(int argc, char *argv[])
             << runTime.cpuTimeIncrement() << " s" << nl << endl;
 
 
-        faces2.writeObject(streamOpt, true);
+        faces2.writeObject
+        (
+            IOstreamOption(format),
+            true
+        );
 
         Info<< "Written new format faceList in = "
             << runTime.cpuTimeIncrement() << " s" << nl << endl;

applications/test/DLList/Test-DLList.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,11 +67,11 @@ int main(int argc, char *argv[])
 
     for (int i = 0; i<10; i++)
     {
-        myList.push_back(1.3*i);
+        myList.append(1.3*i);
     }
 
-    myList.push_back(100.3);
-    myList.push_back(500.3);
+    myList.append(100.3);
+    myList.append(500.3);
 
     Info<< "DLList<scalar>" << nl;
     Info<< nl << "flat-output: " << flatOutput(myList) << nl;
@@ -144,23 +144,23 @@ int main(int argc, char *argv[])
         Info<< " => " << flatOutput(myList) << nl;
     }
 
-    myList.push_back(500.3);
-    myList.push_back(200.3);
-    myList.push_back(100.3);
+    myList.append(500.3);
+    myList.append(200.3);
+    myList.append(100.3);
 
     Info<< nl << "Testing swapUp and swapDown:" << nl;
     Info<< " => " << flatOutput(myList) << nl;
 
     {
-        myList.swapUp(myList.DLListBase::front());
-        myList.swapUp(myList.DLListBase::back());
+        myList.swapUp(myList.DLListBase::first());
+        myList.swapUp(myList.DLListBase::last());
 
         Info<< nl << "swapUp => " << flatOutput(myList) << nl;
     }
 
     {
-        myList.swapDown(myList.DLListBase::front());
-        myList.swapDown(myList.DLListBase::back());
+        myList.swapDown(myList.DLListBase::first());
+        myList.swapDown(myList.DLListBase::last());
 
         Info<< nl << "swapDown => " << flatOutput(myList) << nl;
     }
@@ -200,7 +200,7 @@ int main(int argc, char *argv[])
 
         for (int i = 0; i<5; i++)
         {
-            labList.push_back(identity(6));
+            labList.append(identity(6));
         }
 
         Info<< nl
@@ -221,16 +221,16 @@ int main(int argc, char *argv[])
         List<label> content1 = identity(10);
 
         Info<< nl
-            << " move push_back ";
+            << " move append ";
         printAddress(content1);
 
-        labList.push_back(std::move(content1));
+        labList.append(std::move(content1));
 
         Info<< " content " << flatOutput(content1) << nl
             << " list" << labList  << nl;
 
         printAddresses(labList);
-        // labList.push_back(content1);
+        // labList.append(content1);
     }
 
     Info<< nl << "Done." << endl;

applications/test/DirLister/DirLister.C

@@ -28,6 +28,11 @@ License
 #include "DirLister.H"
 #include <dirent.h>
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+static const Foam::word extgz("gz");
+
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 bool Foam::DirLister::const_iterator::open(const fileName& dir)
@@ -105,9 +110,9 @@ Foam::word Foam::DirLister::next(DIR* dirPtr) const
 
         if (ok)
         {
-            if (fType == fileName::FILE && stripgz_ && name.has_ext("gz"))
+            if (fType == fileName::FILE && stripgz_ && name.hasExt(extgz))
             {
-                name.remove_ext();
+                name = name.lessExt();
             }
 
             if (!name.empty() && accept(name))

applications/test/DynamicList2/Test-DynamicList2.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -104,7 +104,7 @@ void readList
 {
     OTstream os;
     os << input;
-    ITstream is(os.tokens());
+    ITstream is("input", os.tokens());
 
     is >> output;
 }
@@ -118,7 +118,7 @@ void readList
 {
     OTstream os;
     os << input;
-    ITstream is(os.tokens());
+    ITstream is("input", os.tokens());
 
     is >> output;
 }

applications/test/FieldFields1/Make/files

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

applications/test/FieldFields1/Test-FieldFields1.C

@@ -1,117 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    Test-FieldFields1
-
-\*---------------------------------------------------------------------------*/
-
-#include "symmTensorField.H"
-#include "tensorField.H"
-#include "FieldFields.H"
-#include "Random.H"
-
-using namespace Foam;
-
-
-template<class Cmpt>
-void printFieldField(const FieldField<Field, Cmpt>& ff)
-{
-    forAll(ff, i)
-    {
-        Info<< i << ": " << flatOutput(ff[i]) << nl;
-    }
-    Info<< nl;
-}
-
-
-template<class Type>
-tmp<Field<Type>> randomField(Random& rnd, label dim)
-{
-    auto tfld = tmp<Field<Type>>::New(dim);
-    auto& fld = tfld.ref();
-
-    for (Type& val : fld)
-    {
-        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
-        {
-            setComponent(val, cmpt) = rnd.position<label>(0, 100);
-        }
-    }
-
-    return tfld;
-}
-
-
-template<class Type>
-tmp<Field<Type>> randomField(label dim)
-{
-    Random rnd;
-    return randomField<Type>(rnd, dim);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    // scalarField
-    {
-        Info<< nl << "scalarFieldField" << nl;
-
-        Random rnd;
-
-        FieldField<Field, scalar> sff1(6);
-        forAll(sff1, i)
-        {
-            sff1.set(i, randomField<scalar>(rnd, 8));
-        }
-
-        printFieldField(sff1);
-
-        Info<< nl << "indexing:" << nl;
-
-        {
-            labelPair index;
-            const label range1 = sff1.size()-1;
-            const label range2 = sff1[0].size()-1;
-
-            for (label iter = 0; iter < 10; ++iter)
-            {
-                index.first() = rnd.position<label>(0, range1);
-                index.second() = rnd.position<label>(0, range2);
-
-                Info<< index << " => " << sff1[index] << nl;
-            }
-        }
-    }
-
-    Info<< nl << "End\n" << nl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/test/FieldFields2/Make/files

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

applications/test/FixedList/Test-FixedList.C

@@ -189,12 +189,6 @@ int main(int argc, char *argv[])
             << " hash:" << FixedList<label, 4>::hasher()(list1) << nl
             << " hash:" << Hash<FixedList<label, 4>>()(list1) << nl;
 
-        Info<< "get<0>: " << list1.get<0>() << nl;
-        Info<< "get<1>: " << list1.get<1>() << nl;
-        Info<< "get<2>: " << list1.get<2>() << nl;
-        Info<< "get<3>: " << list1.get<3>() << nl;
-// Will not compile:  Info<< "get<4>: " << list1.get<4>() << nl;
-
         label a[4] = {0, 1, 2, 3};
         FixedList<label, 4> list2(a);
 

applications/test/Function1/case1/constant/function1Properties

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2212                                 |
+|  \\    /   O peration     | Version:  v2206                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/