ENH: regIOobject: factory method from IO.

- compiles - deepClone (not needed at the moment) - objectRegistry::subRegistry with fileName - Time::null object so IOobject can be constructed without objectRegistry
2024-01-11 15:15:34 +01:00
14504 changed files with 157478 additions and 318637 deletions
--- a/.gitlab/issue_templates/bug.md
+++ b/.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 : v1806|v1812|v1906|v1912|v2006|v2012|v2106 etc
  Operating system : ubuntu|openSUSE|centos etc
  Hardware info    : any info that may help?
  Compiler         : gcc|intel|clang etc
--- a/13
+++ b/13
@ -37,7 +37,7 @@ echo
 echo ========================================
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"
 echo "Starting compile ${WM_PROJECT_DIR##*/} ${0##*/}"
-echo "  $WM_COMPILER ${WM_COMPILER_TYPE:-system} compiler [${WM_COMPILE_CONTROL}]"
+echo "  $WM_COMPILER ${WM_COMPILER_TYPE:-system} compiler"
 echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
 echo ========================================
 echo
@ -70,16 +70,15 @@ applications/Allwmake $targetType $*


 # Additional components/modules
-case "$FOAM_MODULE_PREFIX" in
-(false | none)
+if [ "$FOAM_MODULE_PREFIX" = false ] || [ "$FOAM_MODULE_PREFIX" = none ]
+then
    echo ========================================
    echo "OpenFOAM modules disabled (prefix=${FOAM_MODULE_PREFIX})"
    echo
-    ;;
-(*)
+else
    # Use wmake -all instead of Allwmake to allow for overrides
-    ( cd "$WM_PROJECT_DIR/modules" 2>/dev/null && wmake -all )
-esac
+    (cd "$WM_PROJECT_DIR/modules" 2>/dev/null && wmake -all)
+fi

 # Count files in given directory. Ignore "Test-*" binaries.
 _foamCountDirEntries()
--- a/CONTRIBUTORS.md
+++ b/CONTRIBUTORS.md
@ -5,8 +5,6 @@ It is likely incomplete...

 ## Contributors (alphabetical by surname)

- Tetsuo Aoyagi
- Akira Azami
 - William Bainbridge
 - Gabriel Barajas
 - Kutalmis Bercin
@ -21,7 +19,6 @@ It is likely incomplete...
 - Bernhard Gschaider
 - Andrew Heather
 - David Hill
- Yoshiaki Inoue
 - Mattijs Janssens
 - Andrew Jackson
 - Hrvoje Jasak
--- a/META-INFO/api-info
+++ b/META-INFO/api-info
@ -1,2 +1,2 @@
-api=2212
+api=2108
 patch=0
--- a/README.md
+++ b/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-v2106 version:
 ```
-source /installation/path/OpenFOAM-v2212/etc/bashrc
+source /installation/path/OpenFOAM-v2106/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-v2106
+\-- ThirdParty-v2106
 ```
 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-sandbox2106, 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*, `v2106-myCustom`,
     without requiring a renamed ThirdParty. The API value would still
-     be `2212` and the original `ThirdParty-v2212/` would be found.
+     be `2106` and the original `ThirdParty-v2106/` would be found.
 4. PREFIX/ThirdParty-API
   * same as the previous example, but using an unadorned API value.
 5. PREFIX/ThirdParty-common
@ -213,4 +213,4 @@ ThirdParty directory will contain either an `Allwmake` file or a
 - [Governance](http://www.openfoam.com/governance/), [Governance Projects](https://www.openfoam.com/governance/projects)
 - [Contacting OpenCFD](http://www.openfoam.com/contact/)

-Copyright 2016-2022 OpenCFD Ltd
+Copyright 2016-2021 OpenCFD Ltd
--- a/applications/Allwmake
+++ b/applications/Allwmake
@ -21,10 +21,10 @@ echo ========================================

 #------------------------------------------------------------------------------

-wmake -all $targetType utilities
 wmake -all $targetType solvers
+wmake -all $targetType utilities

-# Extra tools (optional)
+# Optional
 ## wmake -all $targetType tools

 #------------------------------------------------------------------------------
--- a/applications/solvers/acoustic/acousticFoam/Make/options
+++ b/applications/solvers/acoustic/acousticFoam/Make/options
@ -11,5 +11,4 @@ EXE_LIBS = \
    -lfvOptions \
    -lmeshTools \
    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/options
+++ b/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 \
--- a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H
+++ b/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)
-        );
-    }
--- a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/overLaplacianDyMFoam.C
+++ b/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"
--- a/applications/solvers/basic/potentialFoam/overPotentialFoam/overPotentialFoam.C
+++ b/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
--- a/applications/solvers/combustion/PDRFoam/Make/options
+++ b/applications/solvers/combustion/PDRFoam/Make/options
@ -25,7 +25,6 @@ EXE_LIBS = \
    -lengine \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
--- a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
@ -116,7 +116,7 @@ int main(int argc, char *argv[])
            "normalisedGradP",
            tmagGradP()/max(tmagGradP())
        );
-        normalisedGradP.writeOpt(IOobject::AUTO_WRITE);
+        normalisedGradP.writeOpt() = IOobject::AUTO_WRITE;
        tmagGradP.clear();

        ++runTime;
--- a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.H
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.H
@ -138,6 +138,13 @@ public:
            const surfaceScalarField& phi
        );

+        //- Return a clone
+        virtual refPtr<regIOobject> deepClone() const
+        {
+            NotImplemented;
+            return nullptr;
+        }
+

    //- Destructor
    virtual ~PDRDragModel();
--- a/applications/solvers/combustion/XiFoam/Make/options
+++ b/applications/solvers/combustion/XiFoam/Make/options
@ -19,7 +19,6 @@ EXE_LIBS = \
    -lengine \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
--- a/applications/solvers/combustion/XiFoam/XiDyMFoam/Make/options
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/Make/options
@ -26,7 +26,6 @@ EXE_LIBS = \
    -lengine \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
--- a/applications/solvers/combustion/XiFoam/XiEngineFoam/Make/options
+++ b/applications/solvers/combustion/XiFoam/XiEngineFoam/Make/options
@ -20,7 +20,6 @@ EXE_LIBS = \
    -lengine \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
--- a/applications/solvers/combustion/coldEngineFoam/Make/options
+++ b/applications/solvers/combustion/coldEngineFoam/Make/options
@ -20,7 +20,6 @@ EXE_LIBS = \
    -lengine \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie
--- a/applications/solvers/combustion/fireFoam/Make/options
+++ b/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 \
@ -34,7 +34,6 @@ EXE_LIBS = \
    -lsampling \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/combustion/fireFoam/createFieldRefs.H
+++ b/applications/solvers/combustion/fireFoam/createFieldRefs.H
@ -1,4 +1,4 @@
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
 regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
-const label inertIndex(composition.species().find(inertSpecie));
+const label inertIndex(composition.species()[inertSpecie]);
--- a/applications/solvers/combustion/reactingFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/Make/options
@ -19,7 +19,6 @@ EXE_LIBS = \
    -lsampling \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lreactionThermophysicalModels \
    -lspecie \
    -lcompressibleTransportModels \
--- a/applications/solvers/combustion/reactingFoam/createFieldRefs.H
+++ b/applications/solvers/combustion/reactingFoam/createFieldRefs.H
@ -1,3 +1,3 @@
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-const label inertIndex(composition.species().find(inertSpecie));
+const label inertIndex(composition.species()[inertSpecie]);
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options
@ -20,7 +20,6 @@ EXE_LIBS = \
    -lsampling \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lreactionThermophysicalModels \
    -lspecie \
    -lcompressibleTransportModels \
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFieldRefs.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFieldRefs.H
@ -1,3 +1,3 @@
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-const label inertIndex(composition.species().find(inertSpecie));
+const label inertIndex(composition.species()[inertSpecie]);
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/Make/options
@ -20,7 +20,6 @@ EXE_LIBS = \
    -lsampling \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lreactionThermophysicalModels \
    -lspecie \
    -lcompressibleTransportModels \
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFieldRefs.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFieldRefs.H
@ -1,3 +1,3 @@
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-const label inertIndex(composition.species().find(inertSpecie));
+const label inertIndex(composition.species()[inertSpecie]);
--- a/applications/solvers/compressible/rhoCentralFoam/Allwclean
+++ b/applications/solvers/compressible/rhoCentralFoam/Allwclean
@ -4,5 +4,6 @@ cd "${0%/*}" || exit    # Run from this directory

 wclean libso BCs
 wclean
+wclean rhoCentralDyMFoam

 #------------------------------------------------------------------------------
--- a/applications/solvers/compressible/rhoCentralFoam/Allwmake
+++ b/applications/solvers/compressible/rhoCentralFoam/Allwmake
@ -5,7 +5,8 @@ cd "${0%/*}" || exit                                # Run from this directory

 (
    wmake $targetType BCs \
- && wmake $targetType
+ && wmake $targetType \
+ && wmake $targetType rhoCentralDyMFoam \
 )

 #------------------------------------------------------------------------------
--- a/applications/solvers/compressible/rhoCentralFoam/Make/options
+++ b/applications/solvers/compressible/rhoCentralFoam/Make/options
@ -7,7 +7,7 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -18,7 +18,4 @@ EXE_LIBS = \
    -lspecie \
    -lrhoCentralFoam \
    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lthermoTools \
-    -ldynamicFvMesh \
-    -ltopoChangerFvMesh
+    -lcompressibleTurbulenceModels
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/files
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/files
@ -0,0 +1,3 @@
+rhoCentralDyMFoam.C
+
+EXE = $(FOAM_APPBIN)/rhoCentralDyMFoam
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options
@ -0,0 +1,26 @@
+EXE_INC = \
+    -I.. \
+    -I../BCs/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lrhoCentralFoam \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
@ -0,0 +1,288 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rhoCentralDyMFoam
+
+Group
+    grpCompressibleSolvers grpMovingMeshSolvers
+
+Description
+    Density-based compressible flow solver based on central-upwind
+    schemes of Kurganov and Tadmor
+    with support for mesh-motion and topology changes.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "fixedRhoFvPatchScalarField.H"
+#include "directionInterpolate.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+#include "motionSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Density-based compressible flow solver based on central-upwind"
+        " schemes of Kurganov and Tadmor.\n"
+        "With support for mesh-motion and topology changes."
+    );
+
+    #define NO_CONTROL
+    #include "postProcess.H"
+
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createTimeControls.H"
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    #include "readFluxScheme.H"
+
+    const dimensionedScalar v_zero(dimVolume/dimTime, Zero);
+
+    // Courant numbers used to adjust the time-step
+    scalar CoNum = 0.0;
+    scalar meanCoNum = 0.0;
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "setDeltaT.H"
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Do any mesh changes
+        mesh.update();
+
+        // --- Directed interpolation of primitive fields onto faces
+
+        surfaceScalarField rho_pos(interpolate(rho, pos));
+        surfaceScalarField rho_neg(interpolate(rho, neg));
+
+        surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
+        surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
+
+        volScalarField rPsi("rPsi", 1.0/psi);
+        surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
+        surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
+
+        surfaceScalarField e_pos(interpolate(e, pos, T.name()));
+        surfaceScalarField e_neg(interpolate(e, neg, T.name()));
+
+        surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
+        surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
+
+        surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
+        surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
+
+        surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
+        surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
+
+        // Make fluxes relative to mesh-motion
+        if (mesh.moving())
+        {
+            phiv_pos -= mesh.phi();
+            phiv_neg -= mesh.phi();
+        }
+        // Note: extracted out the orientation so becomes unoriented
+        phiv_pos.setOriented(false);
+        phiv_neg.setOriented(false);
+
+        volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
+        surfaceScalarField cSf_pos
+        (
+            "cSf_pos",
+            interpolate(c, pos, T.name())*mesh.magSf()
+        );
+        surfaceScalarField cSf_neg
+        (
+            "cSf_neg",
+            interpolate(c, neg, T.name())*mesh.magSf()
+        );
+
+        surfaceScalarField ap
+        (
+            "ap",
+            max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
+        );
+        surfaceScalarField am
+        (
+            "am",
+            min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
+        );
+
+        surfaceScalarField a_pos("a_pos", ap/(ap - am));
+
+        surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
+
+        surfaceScalarField aSf("aSf", am*a_pos);
+
+        if (fluxScheme == "Tadmor")
+        {
+            aSf = -0.5*amaxSf;
+            a_pos = 0.5;
+        }
+
+        surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
+
+        phiv_pos *= a_pos;
+        phiv_neg *= a_neg;
+
+        surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
+        surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
+
+        // Reuse amaxSf for the maximum positive and negative fluxes
+        // estimated by the central scheme
+        amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
+
+        #include "centralCourantNo.H"
+
+        phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
+
+        surfaceVectorField phiU(aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg);
+        // Note: reassembled orientation from the pos and neg parts so becomes
+        // oriented
+        phiU.setOriented(true);
+
+        surfaceVectorField phiUp(phiU + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf());
+
+        surfaceScalarField phiEp
+        (
+            "phiEp",
+            aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
+          + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
+          + aSf*p_pos - aSf*p_neg
+        );
+
+        // Make flux for pressure-work absolute
+        if (mesh.moving())
+        {
+            surfaceScalarField phia(a_pos*p_pos + a_neg*p_neg);
+            phia.setOriented(true);
+
+            phiEp += mesh.phi()*phia;
+        }
+
+        volScalarField muEff("muEff", turbulence->muEff());
+        volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
+
+        // --- Solve density
+        solve(fvm::ddt(rho) + fvc::div(phi));
+
+        // --- Solve momentum
+        solve(fvm::ddt(rhoU) + fvc::div(phiUp));
+
+        U.ref() =
+            rhoU()
+           /rho();
+        U.correctBoundaryConditions();
+        rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();
+
+        if (!inviscid)
+        {
+            solve
+            (
+                fvm::ddt(rho, U) - fvc::ddt(rho, U)
+              - fvm::laplacian(muEff, U)
+              - fvc::div(tauMC)
+            );
+            rhoU = rho*U;
+        }
+
+        // --- Solve energy
+        surfaceScalarField sigmaDotU
+        (
+            "sigmaDotU",
+            (
+                fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
+              + fvc::dotInterpolate(mesh.Sf(), tauMC)
+            )
+          & (a_pos*U_pos + a_neg*U_neg)
+        );
+
+        solve
+        (
+            fvm::ddt(rhoE)
+          + fvc::div(phiEp)
+          - fvc::div(sigmaDotU)
+        );
+
+        e = rhoE/rho - 0.5*magSqr(U);
+        e.correctBoundaryConditions();
+        thermo.correct();
+        rhoE.boundaryFieldRef() ==
+            rho.boundaryField()*
+            (
+                e.boundaryField() + 0.5*magSqr(U.boundaryField())
+            );
+
+        if (!inviscid)
+        {
+            solve
+            (
+                fvm::ddt(rho, e) - fvc::ddt(rho, e)
+              - fvm::laplacian(turbulence->alphaEff(), e)
+            );
+            thermo.correct();
+            rhoE = rho*(e + 0.5*magSqr(U));
+        }
+
+        p.ref() =
+            rho()
+           /psi();
+        p.correctBoundaryConditions();
+        rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
+
+        turbulence->correct();
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -31,14 +30,12 @@ Group
    grpCompressibleSolvers

 Description
-    Density-based compressible flow solver based on
-    central-upwind schemes of Kurganov and Tadmor with
-    support for mesh-motion and topology changes.
+    Density-based compressible flow solver based on central-upwind
+    schemes of Kurganov and Tadmor.

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

 #include "fvCFD.H"
-#include "dynamicFvMesh.H"
 #include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "fixedRhoFvPatchScalarField.H"
@ -52,9 +49,8 @@ int main(int argc, char *argv[])
 {
    argList::addNote
    (
-        "Density-based compressible flow solver based on"
-        " central-upwind schemes of Kurganov and Tadmor with"
-        " support for mesh-motion and topology changes."
+        "Density-based compressible flow solver based on central-upwind"
+        " schemes of Kurganov and Tadmor."
    );

    #define NO_CONTROL
@ -63,7 +59,7 @@ int main(int argc, char *argv[])
    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
+    #include "createMesh.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
    #include "createTimeControls.H"
@ -84,18 +80,6 @@ int main(int argc, char *argv[])

    while (runTime.run())
    {
-        #include "readTimeControls.H"
-
-        if (!LTS)
-        {
-            #include "setDeltaT.H"
-
-            ++runTime;
-
-            // Do any mesh changes
-            mesh.update();
-        }
-
        // --- Directed interpolation of primitive fields onto faces

        surfaceScalarField rho_pos(interpolate(rho, pos));
@ -123,15 +107,6 @@ int main(int argc, char *argv[])
        surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
        phiv_neg.setOriented(false);

-        // Make fluxes relative to mesh-motion
-        if (mesh.moving())
-        {
-            surfaceScalarField meshPhi(mesh.phi());
-            meshPhi.setOriented(false);
-            phiv_pos -= meshPhi;
-            phiv_neg -= meshPhi;
-        }
-
        volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
        surfaceScalarField cSf_pos
        (
@ -182,13 +157,18 @@ int main(int argc, char *argv[])
        amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));

        #include "centralCourantNo.H"
+        #include "readTimeControls.H"

        if (LTS)
        {
            #include "setRDeltaT.H"
-
-            ++runTime;
        }
+        else
+        {
+            #include "setDeltaT.H"
+        }
+
+        ++runTime;

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

@ -209,14 +189,6 @@ int main(int argc, char *argv[])
          + aSf*p_pos - aSf*p_neg
        );

-        // Make flux for pressure-work absolute
-        if (mesh.moving())
-        {
-            surfaceScalarField meshPhi(mesh.phi());
-            meshPhi.setOriented(false);
-            phiEp += meshPhi*(a_pos*p_pos + a_neg*p_neg);
-        }
-
        volScalarField muEff("muEff", turbulence->muEff());
        volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));

--- a/applications/solvers/compressible/rhoPimpleAdiabaticFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleAdiabaticFoam/Make/options
@ -6,7 +6,7 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \

 EXE_LIBS = \
    -lfiniteVolume \
@ -17,5 +17,4 @@ EXE_LIBS = \
    -lfluidThermophysicalModels \
    -lspecie \
    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lthermoTools
+    -lcompressibleTurbulenceModels
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@ -20,11 +20,9 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -lsampling \
    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/Make/options
@ -21,7 +21,6 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -loverset
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createControls.H
@ -0,0 +1,4 @@
+bool ddtCorr
+(
+    pimple.dict().getOrDefault("ddtCorr", true)
+);
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/createFields.H
+++ b/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
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/overRhoPimpleDyMFoam.C
+++ b/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"
                }
            }

--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
+++ b/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;
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/readControls.H
+++ b/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);
--- a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
@ -49,7 +49,7 @@
        );
    }

-    // Update the boundary values of the reciprocal time-step
+    // Update tho boundary values of the reciprocal time-step
    rDeltaT.correctBoundaryConditions();

    Info<< "Flow time scale min/max = "
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@ -19,7 +19,5 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/options
@ -21,7 +21,6 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -lsampling \
--- a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H
@ -1,52 +1,22 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
-------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+    Info<< "Create dynamic mesh for time = "
+        << runTime.timeName() << nl << endl;

-Description
-    Create a dynamicFvMesh, with init()
-
-Required Variables
-    - runTime [Time]
-
-Provided Variables
-    - mesh [dynamicFvMesh], meshPtr
-
-\*---------------------------------------------------------------------------*/
-
-Foam::autoPtr<Foam::dynamicFvMesh> meshPtr;
-
-{
-    Foam::Info << "Create dynamic mesh for time = "
-        << runTime.timeName() << Foam::nl;
-
-    meshPtr = dynamicFvMesh::New
+    autoPtr<dynamicFvMesh> meshPtr
    (
-        IOobject
+        dynamicFvMesh::New
        (
-            polyMesh::defaultRegion,
-            runTime.timeName(),
-            runTime,
-            IOobject::MUST_READ
+            IOobject
+            (
+                polyMesh::defaultRegion,
+                runTime.timeName(),
+                runTime,
+                IOobject::MUST_READ
+            )
        )
    );
-}

+    dynamicFvMesh& mesh = meshPtr();

-dynamicFvMesh& mesh = meshPtr();
-
-// Calculate initial mesh-to-mesh mapping. Note that this should be
-// done under the hood, e.g. as a MeshObject
-mesh.update();
-
-Foam::Info << Foam::endl;
-
-
-// ************************************************************************* //
+    // Calculate initial mesh-to-mesh mapping. Note that this should be
+    // done under the hood, e.g. as a MeshObject
+    mesh.update();
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/Make/options
@ -18,5 +18,4 @@ EXE_LIBS = \
    -lfluidThermophysicalModels \
    -lspecie \
    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lthermoTools
+    -lcompressibleTurbulenceModels
--- a/applications/solvers/compressible/sonicFoam/Make/options
+++ b/applications/solvers/compressible/sonicFoam/Make/options
@ -16,5 +16,4 @@ EXE_LIBS = \
    -lfluidThermophysicalModels \
    -lspecie \
    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lthermoTools
+    -lcompressibleTurbulenceModels
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options
@ -21,7 +21,6 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
--- a/applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H
+++ b/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();
--- a/applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H
+++ b/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);
--- a/applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H
+++ b/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
    (
--- a/applications/solvers/finiteArea/surfactantFoam/createVolFields.H
+++ b/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
    (
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
@ -1,11 +1,8 @@
 EXE_INC = \
    -I../buoyantBoussinesqSimpleFoam \
-    -I../../incompressible/pimpleFoam \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
@ -16,8 +13,6 @@ EXE_LIBS = \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
-    -ldynamicMesh \
-    -ldynamicFvMesh \
    -lsampling \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -31,8 +30,7 @@ Group
    grpHeatTransferSolvers

 Description
-    Transient solver for buoyant, turbulent flow of incompressible fluids,
-    with optional mesh motion and mesh topology changes.
+    Transient solver for buoyant, turbulent flow of incompressible fluids.

    Uses the Boussinesq approximation:
    \f[
@ -53,11 +51,9 @@ Description
 \*---------------------------------------------------------------------------*/

 #include "fvCFD.H"
-#include "dynamicFvMesh.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "radiationModel.H"
-#include "CorrectPhi.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"

@ -68,8 +64,7 @@ int main(int argc, char *argv[])
    argList::addNote
    (
        "Transient solver for buoyant, turbulent flow"
-        " of incompressible fluids, with optional mesh"
-        " motion and mesh topology changes.\n"
+        " of incompressible fluids.\n"
        "Uses the Boussinesq approximation."
    );

@ -78,10 +73,10 @@ int main(int argc, char *argv[])
    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
-    #include "createDyMControls.H"
+    #include "createMesh.H"
+    #include "createControl.H"
    #include "createFields.H"
-    #include "createUfIfPresent.H"
+    #include "createTimeControls.H"
    #include "CourantNo.H"
    #include "setInitialDeltaT.H"
    #include "initContinuityErrs.H"
@ -94,7 +89,7 @@ int main(int argc, char *argv[])

    while (runTime.run())
    {
-        #include "readDyMControls.H"
+        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"

@ -105,34 +100,6 @@ int main(int argc, char *argv[])
        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
-            if (pimple.firstIter() || moveMeshOuterCorrectors)
-            {
-                // Do any mesh changes
-                mesh.controlledUpdate();
-
-                if (mesh.changing())
-                {
-                    MRF.update();
-
-                    if (correctPhi)
-                    {
-                        // Calculate absolute flux
-                        // from the mapped surface velocity
-                        phi = mesh.Sf() & Uf();
-
-                        #include "correctPhi.H"
-
-                        // Make the flux relative to the mesh motion
-                        fvc::makeRelative(phi, U);
-                    }
-
-                    if (checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
-                }
-            }
-
            #include "UEqn.H"
            #include "TEqn.H"

--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@ -42,12 +42,6 @@
            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
-
-            // Correct Uf if the mesh is moving
-            fvc::correctUf(Uf, U, phi);
-
-            // Make the fluxes relative to the mesh motion
-            fvc::makeRelative(phi, U);
        }
    }

--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../../compressible/rhoPimpleFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
@ -8,8 +7,6 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/regionFaModels\lnInclude

 EXE_LIBS = \
@ -23,10 +20,5 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -latmosphericModels \
-    -ldynamicMesh \
-    -ldynamicFvMesh \
-    -ltopoChangerFvMesh \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -31,9 +30,8 @@ Group
    grpHeatTransferSolvers

 Description
-    Transient solver for buoyant, turbulent flow of compressible fluids
-    for ventilation and heat-transfer, with optional mesh motion
-    and mesh topology changes.
+    Transient solver for buoyant, turbulent flow of compressible fluids for
+    ventilation and heat-transfer.

    Turbulence is modelled using a run-time selectable compressible RAS or
    LES model.
@ -41,16 +39,12 @@ Description
 \*---------------------------------------------------------------------------*/

 #include "fvCFD.H"
-#include "dynamicFvMesh.H"
 #include "rhoThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "radiationModel.H"
-#include "CorrectPhi.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
 #include "pressureControl.H"
-#include "localEulerDdtScheme.H"
-#include "fvcSmooth.H"

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

@ -59,8 +53,7 @@ int main(int argc, char *argv[])
    argList::addNote
    (
        "Transient solver for buoyant, turbulent fluid flow"
-        " of compressible fluids, including radiation,"
-        " with optional mesh motion and mesh topology changes."
+        " of compressible fluids, including radiation."
    );

    #include "postProcess.H"
@ -68,105 +61,36 @@ int main(int argc, char *argv[])
    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
-    #include "createDyMControls.H"
-    #include "initContinuityErrs.H"
+    #include "createMesh.H"
+    #include "createControl.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
-    #include "createRhoUfIfPresent.H"
+    #include "initContinuityErrs.H"
+    #include "createTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.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 (correctPhi)
-        {
-            divrhoU.reset
-            (
-                new volScalarField
-                (
-                    "divrhoU",
-                    fvc::div(fvc::absolute(phi, rho, U))
-                )
-            );
-        }
-
-        if (LTS)
-        {
-            #include "setRDeltaT.H"
-        }
-        else
-        {
-            #include "compressibleCourantNo.H"
-            #include "setDeltaT.H"
-        }
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"

        ++runTime;

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

+        #include "rhoEqn.H"
+
        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
-            if (pimple.firstIter() || moveMeshOuterCorrectors)
-            {
-                // Store momentum to set rhoUf for introduced faces.
-                autoPtr<volVectorField> rhoU;
-                if (rhoUf.valid())
-                {
-                    rhoU.reset(new volVectorField("rhoU", rho*U));
-                }
-
-                // Do any mesh changes
-                mesh.update();
-
-                if (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 "correctPhi.H"
-
-                        // Make the fluxes relative to the mesh-motion
-                        fvc::makeRelative(phi, rho, U);
-                    }
-
-                    if (checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
-                }
-            }
-
-            if (pimple.firstIter() && !pimple.SIMPLErho())
-            {
-                #include "rhoEqn.H"
-            }
-
            #include "UEqn.H"
            #include "EEqn.H"

--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@ -1,5 +1,3 @@
-#include "createRDeltaT.H"
-
 Info<< "Reading thermophysical properties\n" << endl;

 autoPtr<rhoThermo> pThermo(rhoThermo::New(mesh));
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/Make/options
@ -24,7 +24,6 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -latmosphericModels \
    -loverset \
    -ldynamicFvMesh \
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H
+++ b/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"
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/overBuoyantPimpleDyMFoam.C
+++ b/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);
                }
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H
+++ b/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;
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/readControls.H
+++ b/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);
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@ -28,8 +28,6 @@ MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 // Update the pressure BCs to ensure flux consistency
 constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);

-fvc::makeRelative(phiHbyA, rho, U);
-
 fvScalarMatrix p_rghDDtEqn
 (
    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
@ -106,15 +104,7 @@ else

 rho = thermo.rho();

-// Correct rhoUf if the mesh is moving
-fvc::correctRhoUf(rhoUf, rho, U, phi);
-
 if (thermo.dpdt())
 {
    dpdt = fvc::ddt(p);
-
-    if (mesh.moving())
-    {
-        dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
-    }
 }
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
@ -20,7 +20,5 @@ EXE_LIBS = \
    -lradiationModels \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@ -37,8 +37,6 @@ EXE_LIBS = \
    -lcombustionModels \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lradiationModels \
    -lregionModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2019,2022 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -70,7 +70,6 @@ int main(int argc, char *argv[])
    #define CREATE_MESH createMeshesPostProcess.H
    #include "postProcess.H"

-    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMeshes.H"
@ -113,19 +112,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 +132,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 +165,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"
                    }

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
@ -30,11 +30,10 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lradiationModels \
    -lfvOptions \
    -lfaOptions \
    -lregionModels \
    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
+
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2019,2022 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -62,7 +62,6 @@ int main(int argc, char *argv[])
    #define CREATE_MESH createMeshesPostProcess.H
    #include "postProcess.H"

-    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMeshes.H"
@ -76,21 +75,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 +98,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 +120,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"
                }

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidMeshes.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidMeshes.H
@ -4,61 +4,21 @@

    forAll(fluidNames, i)
    {
-        if (args.dryRun() || args.found("dry-run-write"))
-        {
-            Info
-                << "Operating in 'dry-run' mode:"
-                << " case will run for 1 time step.  "
-                << "All checks assumed OK on a clean exit" << endl;
+        Info<< "Create fluid mesh for region " << fluidNames[i]
+            << " for time = " << runTime.timeName() << nl << endl;

-            FieldBase::allowConstructFromLargerSize = true;
-
-            // Create a simplified 1D mesh and attempt to re-create boundary
-            // conditions
-            fluidRegions.set
+        fluidRegions.set
+        (
+            i,
+            new fvMesh
            (
-                i,
-                new simplifiedMeshes::columnFvMesh(runTime, fluidNames[i])
-            );
-
-            // Stop after 1 iteration of the simplified mesh
-
-            if (args.found("dry-run-write"))
-            {
-                // Using saWriteNow triggers function objects execute(), write()
-                runTime.stopAt(Time::saWriteNow);
-
-                // Make sure mesh gets output to the current time (since
-                // instance no longer constant)
-                fluidRegions[i].setInstance(runTime.timeName());
-            }
-            else
-            {
-                // Using saNoWriteNow triggers function objects execute(),
-                // but not write()
-                runTime.stopAt(Time::saNoWriteNow);
-            }
-
-            functionObject::outputPrefix = "postProcessing-dry-run";
-        }
-        else
-        {
-            Info<< "Create fluid mesh for region " << fluidNames[i]
-                << " for time = " << runTime.timeName() << nl << endl;
-
-            fluidRegions.set
-            (
-                i,
-                new fvMesh
+                IOobject
                (
-                    IOobject
-                    (
-                        fluidNames[i],
-                        runTime.timeName(),
-                        runTime,
-                        IOobject::MUST_READ
-                    )
+                    fluidNames[i],
+                    runTime.timeName(),
+                    runTime,
+                    IOobject::MUST_READ
                )
-            );
-        }
+            )
+        );
    }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H
@ -5,5 +5,3 @@

    const bool momentumPredictor =
        simple.getOrDefault("momentumPredictor", true);
-
-    simple.readIfPresent("frozenFlow", frozenFlowFluid[i]);
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/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");

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/chtMultiRegionTwoPhaseEulerFoam.C
+++ b/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.
@ -67,7 +67,6 @@ int main(int argc, char *argv[])
    #define CREATE_MESH createMeshesPostProcess.H
    #include "postProcess.H"

-    #include "addCheckCaseOptions.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMeshes.H"
@ -108,23 +107,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 +134,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"
                    }
                }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/readFluidMultiRegionPIMPLEControls.H
@ -9,5 +9,3 @@
    (
        pimpleDict.getOrDefault<int>("nEnergyCorrectors", 1)
    );
-
-    pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H
@ -1,3 +1,5 @@
+    fvMesh& mesh = fluidRegions[i];
+
    twoPhaseSystem& fluid  = phaseSystemFluid[i];

    phaseModel& phase1 = fluid.phase1();
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H
@ -4,61 +4,21 @@

    forAll(fluidNames, i)
    {
-        if (args.dryRun() || args.found("dry-run-write"))
-        {
-            Info
-                << "Operating in 'dry-run' mode:"
-                << " case will run for 1 time step.  "
-                << "All checks assumed OK on a clean exit" << endl;
+        Info<< "Create fluid mesh for region " << fluidNames[i]
+            << " for time = " << runTime.timeName() << nl << endl;

-            FieldBase::allowConstructFromLargerSize = true;
-
-            // Create a simplified 1D mesh and attempt to re-create boundary
-            // conditions
-            fluidRegions.set
+        fluidRegions.set
+        (
+            i,
+            new fvMesh
            (
-                i,
-                new simplifiedMeshes::columnFvMesh(runTime, fluidNames[i])
-            );
-
-            // Stop after 1 iteration of the simplified mesh
-
-            if (args.found("dry-run-write"))
-            {
-                // Using saWriteNow triggers function objects execute(), write()
-                runTime.stopAt(Time::saWriteNow);
-
-                // Make sure mesh gets output to the current time (since
-                // instance no longer constant)
-                fluidRegions[i].setInstance(runTime.timeName());
-            }
-            else
-            {
-                // Using saNoWriteNow triggers function objects execute(),
-                // but not write()
-                runTime.stopAt(Time::saNoWriteNow);
-            }
-
-            functionObject::outputPrefix = "postProcessing-dry-run";
-        }
-        else
-        {
-            Info<< "Create fluid mesh for region " << fluidNames[i]
-                << " for time = " << runTime.timeName() << nl << endl;
-
-            fluidRegions.set
-            (
-                i,
-                new fvMesh
+                IOobject
                (
-                    IOobject
-                    (
-                        fluidNames[i],
-                        runTime.timeName(),
-                        runTime,
-                        IOobject::MUST_READ
-                    )
+                    fluidNames[i],
+                    runTime.timeName(),
+                    runTime,
+                    IOobject::MUST_READ
                )
-            );
-        }
+            )
+        );
    }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H
@ -8,5 +8,3 @@

    const bool momentumPredictor =
        pimple.getOrDefault("momentumPredictor", true);
-
-    pimple.readIfPresent("frozenFlow", frozenFlowFluid[i]);
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@ -1,3 +1,5 @@
+    fvMesh& mesh = fluidRegions[i];
+
    CombustionModel<rhoReactionThermo>& reaction = reactionFluid[i];

    rhoReactionThermo& thermo = reaction.thermo();
@ -10,15 +12,15 @@
    if (Y.size())
    {
        const word inertSpecie(thermo.get<word>("inertSpecie"));
-        inertIndex = composition.species().find(inertSpecie);
-        if (inertIndex < 0)
+        if (!composition.species().found(inertSpecie))
        {
            FatalIOErrorIn(args.executable().c_str(), thermo)
                << "Inert specie " << inertSpecie
                << " not found in available species "
-                << flatOutput(composition.species())
+                << composition.species()
                << exit(FatalIOError);
        }
+        inertIndex = composition.species()[inertSpecie];
    }

    volScalarField& rho = rhoFluid[i];
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
@ -35,7 +35,7 @@
                (
                    solidRegions[i],
                    thermos[i],
-                    coordinateSystem::typeName
+                    coordinateSystem::typeName_()
                )
            );

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H
@ -4,67 +4,26 @@

    forAll(solidNames, i)
    {
-        if (args.dryRun() || args.found("dry-run-write"))
-        {
-            Info
-                << "Operating in 'dry-run' mode:"
-                << " case will run for 1 time step.  "
-                << "All checks assumed OK on a clean exit" << endl;
+        Info<< "Create solid mesh for region " << solidNames[i]
+            << " for time = " << runTime.timeName() << nl << endl;

-            FieldBase::allowConstructFromLargerSize = true;
-
-            // Create a simplified 1D mesh and attempt to re-create boundary
-            // conditions
-            solidRegions.set
+        solidRegions.set
+        (
+            i,
+            new fvMesh
            (
-                i,
-                new simplifiedMeshes::columnFvMesh(runTime, solidNames[i])
-            );
-
-            // Stop after 1 iteration of the simplified mesh
-
-            if (args.found("dry-run-write"))
-            {
-                // Using saWriteNow triggers function objects execute(), write()
-                runTime.stopAt(Time::saWriteNow);
-
-                // Make sure mesh gets output to the current time (since
-                // instance no longer constant)
-                solidRegions[i].setInstance(runTime.timeName());
-            }
-            else
-            {
-                // Using saNoWriteNow triggers function objects execute(),
-                // but not write()
-                runTime.stopAt(Time::saNoWriteNow);
-            }
-
-            functionObject::outputPrefix = "postProcessing-dry-run";
-        }
-        else
-        {
-            Info<< "Create solid mesh for region " << solidNames[i]
-                << " for time = " << runTime.timeName() << nl << endl;
-
-            solidRegions.set
-            (
-                i,
-                new fvMesh
+                IOobject
                (
-                    IOobject
-                    (
-                        solidNames[i],
-                        runTime.timeName(),
-                        runTime,
-                        IOobject::MUST_READ
-                    )
+                    solidNames[i],
+                    runTime.timeName(),
+                    runTime,
+                    IOobject::MUST_READ
                )
-            );
+            )
+        );

-            // Force calculation of geometric properties to prevent it being
-            // done
-            // later in e.g. some boundary evaluation
-            //(void)solidRegions[i].weights();
-            //(void)solidRegions[i].deltaCoeffs();
-        }
+        // Force calculation of geometric properties to prevent it being done
+        // later in e.g. some boundary evaluation
+        //(void)solidRegions[i].weights();
+        //(void)solidRegions[i].deltaCoeffs();
    }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
@ -1,3 +1,4 @@
+fvMesh& mesh = solidRegions[i];
 solidThermo& thermo = thermos[i];

 tmp<volScalarField> trho = thermo.rho();
--- a/applications/solvers/heatTransfer/solidFoam/Make/files
+++ b/applications/solvers/heatTransfer/solidFoam/Make/files
@ -1,4 +1,3 @@
-solidRegionDiffNo.C
 solidFoam.C

 EXE = $(FOAM_APPBIN)/solidFoam
--- a/applications/solvers/heatTransfer/solidFoam/Make/options
+++ b/applications/solvers/heatTransfer/solidFoam/Make/options
@ -1,6 +1,4 @@
 EXE_INC = \
-    -I../chtMultiRegionFoam/include \
-    -I../chtMultiRegionFoam/solid \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
@ -23,5 +21,4 @@ EXE_LIBS = \
    -lfluidThermophysicalModels \
    -lsolidThermo \
    -lradiationModels \
-    -lspecie \
-    -lthermoTools
+    -lspecie
--- a/applications/solvers/heatTransfer/solidFoam/createFields.H
+++ b/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();
@ -90,24 +90,5 @@ else
    );
 }

-// Consider mesh flux to correct for mesh deformation
-bool meshFluxCorr(false);
-if (mesh.solutionDict().found("SIMPLE"))
-{
-    meshFluxCorr =
-        mesh.solutionDict().subDict("SIMPLE").getOrDefault<bool>
-        (
-            "meshFluxCorrection", false
-        );
-}
-else if (mesh.solutionDict().found("PIMPLE"))
-{
-    meshFluxCorr =
-        mesh.solutionDict().subDict("PIMPLE").getOrDefault<bool>
-        (
-            "meshFluxCorrection", false
-        );
-}
-
 #include "createRadiationModel.H"
 #include "createFvOptions.H"
--- a/applications/solvers/heatTransfer/solidFoam/dummyCourantNo.H
+++ b/applications/solvers/heatTransfer/solidFoam/dummyCourantNo.H
@ -1,36 +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/>.
-
-Global
-    CourantNo
-
-Description
-    Dummy (fluid) Courant number.
-
-\*---------------------------------------------------------------------------*/
-
-scalar CoNum = 0.0;
-
-// ************************************************************************* //
--- a/applications/solvers/heatTransfer/solidFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/solidFoam/hEqn.H
@ -11,14 +11,14 @@
        fvOptions(rho, h)
    );

-    if (meshFluxCorr)
+    if (mesh.changing())
    {
        surfaceScalarField phihMesh
        (
            fvc::interpolate(betav*rho*h)*mesh.phi()
        );

-       hEqn -= fvc::div(phihMesh);
+        hEqn -= fvc::div(phihMesh);
    }

    hEqn.relax();
--- a/applications/solvers/heatTransfer/solidFoam/solidDiffusionNo.H
+++ b/applications/solvers/heatTransfer/solidFoam/solidDiffusionNo.H
@ -1,31 +0,0 @@
-scalar DiNum = -GREAT;
-
-{
-    tmp<volScalarField> magKappa;
-    if (thermo.isotropic())
-    {
-        magKappa = thermo.kappa();
-    }
-    else
-    {
-        magKappa = mag(thermo.Kappa());
-    }
-
-    tmp<volScalarField> tcp = thermo.Cp();
-    const volScalarField& cp = tcp();
-
-    tmp<volScalarField> trho = thermo.rho();
-    const volScalarField& rho = trho();
-
-    DiNum = max
-    (
-        solidRegionDiffNo
-        (
-            mesh,
-            runTime,
-            rho*cp,
-            magKappa()
-        ),
-        DiNum
-    );
-}
--- a/applications/solvers/heatTransfer/solidFoam/solidFoam.C
+++ b/applications/solvers/heatTransfer/solidFoam/solidFoam.C
@ -41,8 +41,6 @@ Description
 #include "fvOptions.H"
 #include "simpleControl.H"
 #include "pimpleControl.H"
-#include "dummyCourantNo.H"
-#include "solidRegionDiffNo.H"
 #include "coordinateSystem.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -90,31 +88,17 @@ int main(int argc, char *argv[])
    {
        pimpleControl pimple(mesh);

-        #include "createDyMControls.H"
-
        while (runTime.run())
        {
-            #include "readDyMControls.H"
-            #include "readSolidTimeControls.H"
-
-            #include "solidDiffusionNo.H"
-            #include "setMultiRegionDeltaT.H"
-
            ++runTime;

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

            while (pimple.loop())
            {
-                if (pimple.firstIter() || moveMeshOuterCorrectors)
+                if (pimple.firstIter())
                {
-                    // Do any mesh changes
-                    mesh.controlledUpdate();
-
-                    if (mesh.changing() && checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
+                    mesh.update();
                }

                while (pimple.correct())
--- a/applications/solvers/heatTransfer/solidFoam/solidRegionDiffNo.C
+++ b/applications/solvers/heatTransfer/solidFoam/solidRegionDiffNo.C
@ -1 +0,0 @@
-../chtMultiRegionFoam/solid/solidRegionDiffNo.C
--- a/applications/solvers/heatTransfer/solidFoam/solidRegionDiffNo.H
+++ b/applications/solvers/heatTransfer/solidFoam/solidRegionDiffNo.H
@ -1 +0,0 @@
-../chtMultiRegionFoam/solid/solidRegionDiffNo.H
--- a/applications/solvers/heatTransfer/thermoFoam/Make/options
+++ b/applications/solvers/heatTransfer/thermoFoam/Make/options
@ -18,5 +18,4 @@ EXE_LIBS = \
    -lradiationModels \
    -lspecie \
    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lthermoTools
+    -lcompressibleTurbulenceModels
--- a/applications/solvers/incompressible/adjointOptimisationFoam/adjointOptimisationFoam.C
+++ b/applications/solvers/incompressible/adjointOptimisationFoam/adjointOptimisationFoam.C
@ -5,8 +5,8 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2021 FOSS GP
+    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2013-2019 FOSS GP
    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@ -71,6 +71,9 @@ int main(int argc, char *argv[])

        // Solve all adjoint equations
        om.solveAdjointEquations();
+
+        // Compute all sensitivities
+        om.computeSensitivities();
    }

    Info<< "End\n" << endl;
--- a/applications/solvers/incompressible/pimpleFoam/Make/options
+++ b/applications/solvers/incompressible/pimpleFoam/Make/options
@ -22,5 +22,4 @@ EXE_LIBS = \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/incompressible/pimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/pimpleFoam/createFields.H
@ -1,5 +1,3 @@
-#include "createRDeltaT.H"
-
 Info<< "Reading field p\n" << endl;
 volScalarField p
 (
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H
+++ b/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());

--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H
@ -1,4 +1,5 @@
 // Solve the Momentum equation
+
 MRF.correctBoundaryVelocity(U);

 tmp<fvVectorMatrix> tUEqn
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H
+++ b/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)
+);
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H
+++ b/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]);
+    }
+}
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/overPimpleDyMFoam.C
+++ b/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
--- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H
+++ b/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;
--- a/Show More
+++ b/Show More
				`@ -1 +0,0 @@`
				`../chtMultiRegionFoam/solid/solidRegionDiffNo.C`