ENH: new SRFrhoSimpleFoam solver

ENH: SRF - refactoring and enabled ramping of rotation
2021-03-16 21:36:43 +00:00 · 2021-03-16 21:36:43 +00:00
21022 changed files with 3179069 additions and 718964 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 : v2312|v2306|v2212|v2206|v2112 etc
+  OpenFOAM version : v1806|v1812|v1906 etc
  Operating system : ubuntu|openSUSE|centos etc
  Hardware info    : any info that may help?
  Compiler         : gcc|intel|clang etc
--- a/.gitmodules
+++ b/.gitmodules
@ -17,6 +17,3 @@
 [submodule "external-solver"]
 	path = modules/external-solver
 	url = https://develop.openfoam.com/Modules/external-solver.git
-[submodule "turbulence-community"]
-	path = modules/turbulence-community
-	url = https://gitlab.com/openfoam/community/tc-turbulence/turbulence-community.git
--- a/52
+++ b/52
@ -8,7 +8,7 @@ wmake -check-dir "$WM_PROJECT_DIR" 2>/dev/null || {
 }
 if [ -f "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments ]
 then  . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments || \
-    echo "Argument parse error"
+    echo "Argument parse error";
 else
    echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
    echo "    Check your OpenFOAM environment and installation"
@ -23,23 +23,23 @@ else
    echo "mpirun=$(command -v mpirun || true)"
 fi
 echo
-# Report compiler information. First non-blank line from --version output
+# Report compiler information
 compiler="$(wmake -show-path-cxx 2>/dev/null || true)"
 if [ -x "$compiler" ]
 then
    echo "compiler=$compiler"
-    "$compiler" --version 2>/dev/null | sed -e '/^$/d;q'
+    "$compiler" --version 2>/dev/null | sed -ne '1p'
 else
    echo "compiler=unknown"
 fi

 echo
-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 "========================================"
 echo

 # Compile tools for wmake
@ -62,24 +62,26 @@ else
    echo "Skip ThirdParty (no directory)"
 fi

-# OpenFOAM libraries
+echo "========================================"
+echo "Compile OpenFOAM libraries"
+echo
 src/Allwmake $targetType $*

-# OpenFOAM applications
+echo "========================================"
+echo "Compile OpenFOAM applications"
+echo
 applications/Allwmake $targetType $*

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

 # Count files in given directory. Ignore "Test-*" binaries.
 _foamCountDirEntries()
@ -91,18 +93,16 @@ _foamCountDirEntries()
 # Some summary information
 echo
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"
-echo ========================================
+echo "========================================"
 echo "  ${WM_PROJECT_DIR##*/}"
 echo "  $WM_COMPILER ${WM_COMPILER_TYPE:-system} compiler"
 echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
 echo
-
-# The api/patch information
-sed -e 's/^/  /; s/=/ = /' ./META-INFO/api-info 2>/dev/null || true
-
-echo "  bin = $(_foamCountDirEntries "$FOAM_APPBIN") entries"
-echo "  lib = $(_foamCountDirEntries "$FOAM_LIBBIN") entries"
+echo "  api   = $(etc/openfoam -show-api 2>/dev/null)"
+echo "  patch = $(etc/openfoam -show-patch 2>/dev/null)"
+echo "  bin   = $(_foamCountDirEntries "$FOAM_APPBIN") entries"
+echo "  lib   = $(_foamCountDirEntries "$FOAM_LIBBIN") entries"
 echo
-echo ========================================
+echo "========================================"

 #------------------------------------------------------------------------------
--- a/CONTRIBUTORS.md
+++ b/CONTRIBUTORS.md
@ -1,40 +1,29 @@
 # Contributors to OpenFOAM

-The following is a list of known contributors to OpenFOAM.
+The following is an list of known contributors to OpenFOAM.
 It is likely incomplete...

 ## Contributors (alphabetical by surname)

- Horacio Aguerre
- Yu Ankun
- Tetsuo Aoyagi
- Akira Azami
- Pete Bachant
 - William Bainbridge
 - Gabriel Barajas
 - Kutalmis Bercin
- Julius Bergmann
 - Ivor Clifford
 - Greg Collecutt
 - Jonathan Cranford
- Santiago Marquez Damian
 - Sergio Ferraris
 - Matej Forman
 - Marian Fuchs
- Gabriel Gerlero
 - Pawan Ghildiyal
 - Chris Greenshields
 - Bernhard Gschaider
 - Andrew Heather
 - David Hill
- Yoshiaki Inoue
 - Mattijs Janssens
 - Andrew Jackson
 - Hrvoje Jasak
 - Alexander Kabat vel Job
 - Thilo Knacke
- Shannon Leakey
- Sergey Lesnik
 - Tommaso Lucchini
 - Graham Macpherson
 - Alexey Matveichev
@ -47,7 +36,6 @@ It is likely incomplete...
 - Victor Olesen
 - Evangelos Papoutsis-Kiachagias
 - Juho Peltola
- Josep Pocurull
 - Johan Roenby
 - Henrik Rusche
 - Bruno Santos
@ -58,15 +46,11 @@ It is likely incomplete...
 - Gavin Tabor
 - Zeljko Tukovic
 - Eugene De Villiers
- Louis Vittoz
- Vuko Vukcevic
 - Yi Wang
 - Norbert Weber
- Gregor Weiss
- Volker Weissmann
 - Henry Weller
 - Niklas Wikstrom
- Guanyang Xue
 - Thorsten Zirwes

+
 <!----------------------------------------------------------------------------->
--- a/META-INFO/README.md
+++ b/META-INFO/README.md
@ -3,9 +3,8 @@
 Meta-information is generally for OpenFOAM internal use only.

 The format, content and meaning may be changed at anytime without
-prior notice.
-If any of these are changed, these are some of places that will need
-to be updated accordingly:
+notice. If any of these are changed, these are some of places that
+will need to be updated accordingly:

  - bin/foamEtcFile
  - bin/tools/foamConfigurePaths
@ -39,7 +38,7 @@ Format: `date +%y%m`

 - 4-digit year-month (YYMM) integer corresponding to the major
  release or in unusual cases an intermediate release.
-  Example, `2106` for the June-2021 release.
+  Example, `1712` for the Dec-2017 release.


 #### patch
@ -79,9 +78,9 @@ the saved information needs synchronization. The command

 ### Notes

-The saved information is split into two separate files.
-The `api-info` file contains more permanent information,
-whereas the `build-info` is more transient in nature.
+The saved information is split into two separate files. The `api-info`
+contains more permanent information, whereas the `build-info` is more
+transient in nature.

 ----
-2021-06-09
+2020-06-23
--- a/META-INFO/api-info
+++ b/META-INFO/api-info
@ -1,2 +1,2 @@
-api=2402
-patch=240220
+api=2012
+patch=0
--- a/README.md
+++ b/README.md
@ -11,7 +11,7 @@ individual and group contributors, integrations
 [governance guided activities](https://www.openfoam.com/governance/).


-## License
+## Copyright

 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
@ -27,11 +27,10 @@ may redistribute files.
 OpenCFD Ltd grants use of its OpenFOAM trademark by Third Parties on a
 licence basis. ESI Group and OpenFOAM Foundation Ltd are currently
 permitted to use the Name and agreed Domain Name. For information on
-trademark use, please refer to the
-[trademark policy guidelines][link trademark].
+trademark use, please refer to the [trademark policy guidelines](http://www.openfoam.com/legal/trademark-policy.php).

 Please [contact OpenCFD](http://www.openfoam.com/contact) if you have
-any questions about the use of the OpenFOAM trademark.
+any questions on the use of the OpenFOAM trademark.

 Violations of the Trademark are monitored, and will be duly prosecuted.

@ -40,9 +39,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-v2312 version:
+For example, for the OpenFOAM-v1912 version:
 ```
-source /installation/path/OpenFOAM-v2312/etc/bashrc
+source /installation/path/OpenFOAM-v1912/etc/bashrc
 ```

 ## Compiling OpenFOAM
@ -56,11 +55,6 @@ guides:
 | [ThirdParty][repo third] | [readme][link third-readme] | [system requirements][link third-require] | [build][link third-build] |


-If you need to modify the versions or locations of ThirdParty
-software, please read how the
-[OpenFOAM configuration][wiki-config] is structured.
-
-
 ## How do I know which version I am currently using?

 The value of the `$WM_PROJECT_DIR` or even `$WM_PROJECT_VERSION` are
@ -69,19 +63,19 @@ not guaranteed to have any correspondence to the OpenFOAM release
 information is embedded into each application. For example, as
 displayed from `blockMesh -help`:
 ```
-Using: OpenFOAM-com (2012) - visit www.openfoam.com
-Build: b830beb5ea-20210429 (patch=210414)
+Using: OpenFOAM-v1812.local (1812) - visit www.openfoam.com
+Build: 65d6551ff7-20190530 (patch=190531)
 Arch:  LSB;label=32;scalar=64
 ```
 This output contains all of the more interesting information that we need:

 | item                  | value         |
 |-----------------------|---------------|
-| version               | com  (eg, local development branch) |
-| api                   | 2012          |
-| commit                | b830beb5ea    |
-| author date           | 20210429      |
-| patch-level           | (20)210414    |
+| version               | v1812.local   |
+| api                   | 1812          |
+| commit                | 65d6551ff7    |
+| author date           | 20190530      |
+| patch-level           | (20)190531    |
 | label/scalar size     | 32/64 bits    |

 The Arch information may also include the `solveScalar` size
@ -96,19 +90,19 @@ when OpenFOAM has not yet been compiled:
 ```
 $ wmake -build-info
 make
-    api = 2012
-    patch = 210414
+    api = 1812
+    patch = 190531
    branch = master
-    build = 308af39136-20210426
+    build = 65d6551ff7-20190530
 ```
 Similar information is available with `foamEtcFile`, using the
 `-show-api` or `-show-patch` options. For example,
 ```
 $ foamEtcFile -show-api
-2012
+1812

 $ foamEtcFile -show-patch
-210414
+190531
 ```
 This output will generally be the easiest to parse for scripts.
 The `$FOAM_API` convenience environment variable may not reflect the
@ -127,8 +121,8 @@ These 3rd-party sources are normally located in a directory parallel
 to the OpenFOAM directory. For example,
 ```
 /path/parent
-|-- OpenFOAM-v2312
-\-- ThirdParty-v2312
+|-- OpenFOAM-v1912
+\-- ThirdParty-v1912
 ```
 There are, however, many cases where this simple convention is inadequate:

@ -136,7 +130,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-sandbox2312, etc..
+  directory name, e.g. openfoam-sandbox1912, 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,14 +150,17 @@ 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*, `v2312-myCustom`,
+   * allows for an updated value of VERSION, *eg*, `v1912-myCustom`,
     without requiring a renamed ThirdParty. The API value would still
-     be `2312` and the original `ThirdParty-v2312/` would be found.
+     be `1912` and the original `ThirdParty-v1912/` would be found.
 4. PREFIX/ThirdParty-API
-   * same as the previous example, but using an unadorned API value.
+   * this is the same as the previous example, but using an unadorned
+     API value. This also makes sense if the chosen version name also
+     uses the unadorned API value in its naming, *eg*,
+     `1912-patch190131`, `1912.19W03`
 5. PREFIX/ThirdParty-common
-   * permits maximum reuse for various versions, for experienced
-     users who are aware of potential version incompatibilities
+   * permits maximum reuse for various versions, for experienced user
+     who are aware of potential version incompatibilities

 If none of these directories are found to be suitable, it reverts to
 using PROJECT/ThirdParty as a dummy location (even if the directory
@ -187,21 +184,18 @@ ThirdParty directory will contain either an `Allwmake` file or a

 <!-- OpenFOAM -->

-[link trademark]: https://www.openfoam.com/opencfd-limited-trade-mark-policy
-
 [repo openfoam]: https://develop.openfoam.com/Development/openfoam/
 [repo third]: https://develop.openfoam.com/Development/ThirdParty-common/

 [link openfoam-readme]: https://develop.openfoam.com/Development/openfoam/blob/develop/README.md
 [link openfoam-issues]: https://develop.openfoam.com/Development/openfoam/blob/develop/doc/BuildIssues.md
+[link openfoam-config]: https://develop.openfoam.com/Development/openfoam/blob/develop/doc/Config.md
 [link openfoam-build]: https://develop.openfoam.com/Development/openfoam/blob/develop/doc/Build.md
 [link openfoam-require]: https://develop.openfoam.com/Development/openfoam/blob/develop/doc/Requirements.md
 [link third-readme]: https://develop.openfoam.com/Development/ThirdParty-common/blob/develop/README.md
 [link third-build]: https://develop.openfoam.com/Development/ThirdParty-common/blob/develop/BUILD.md
 [link third-require]: https://develop.openfoam.com/Development/ThirdParty-common/blob/develop/Requirements.md

-[wiki-config]: https://develop.openfoam.com/Development/openfoam/-/wikis/configuring
-

 ## Useful Links

@ -209,8 +203,8 @@ ThirdParty directory will contain either an `Allwmake` file or a
 - [Documentation](http://www.openfoam.com/documentation)
 - [Reporting bugs/issues/feature requests](http://www.openfoam.com/code/bug-reporting.php)
 - [Issue tracker](https://develop.openfoam.com/Development/openfoam/-/issues)
- [Code wiki](https://develop.openfoam.com/Development/openfoam/-/wikis/) and [General wiki](http://wiki.openfoam.com/)
- [Governance](http://www.openfoam.com/governance/), [Governance Projects](https://www.openfoam.com/governance/projects)
+- [Code wiki](https://develop.openfoam.com/Development/openfoam/-/wikis/) and [general wiki](http://wiki.openfoam.com/)
+- [Community](http://www.openfoam.com/community/), [Governance](http://www.openfoam.com/governance/)
 - [Contacting OpenCFD](http://www.openfoam.com/contact/)

-Copyright 2016-2023 OpenCFD Ltd
+Copyright 2016-2020 OpenCFD Ltd
--- a/applications/Allwmake
+++ b/applications/Allwmake
@ -8,29 +8,19 @@ wmake -check-dir "$WM_PROJECT_DIR/applications" 2>/dev/null || {
 }
 if [ -f "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments ]
 then  . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments || \
-    echo "Argument parse error"
+    echo "Argument parse error";
 else
    echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
    echo "    Check your OpenFOAM environment and installation"
    exit 1
 fi

-echo ========================================
-echo Compile OpenFOAM applications
-echo ========================================
-
 #------------------------------------------------------------------------------

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

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

 #------------------------------------------------------------------------------
-
-echo ========================================
-echo Done OpenFOAM applications
-echo ========================================
-
-#------------------------------------------------------------------------------
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@ -118,7 +118,7 @@ int main(int argc, char *argv[])
                fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
            );

-            pEqn.solve(p.select(piso.finalInnerIter()));
+            pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));

            phi = phiHbyA - pEqn.flux();

--- 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/acoustic/acousticFoam/readTransportProperties.H
+++ b/applications/solvers/acoustic/acousticFoam/readTransportProperties.H
@ -18,6 +18,6 @@ dimensionedScalar rho("rho", dimDensity, transportProperties);

 scalar MaxCo =
    max(mesh.surfaceInterpolation::deltaCoeffs()*c0).value()
-   *runTime.deltaTValue();
+   *runTime.deltaT().value();

 Info<< "Max acoustic Courant Number = " << MaxCo << endl;
--- 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/createFields.H
+++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H
@ -1,4 +1,125 @@
-#include "../createFields.H"
+Info<< "Reading velocity field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Initialise the velocity internal field to zero
+U = dimensionedVector(U.dimensions(), Zero);
+
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    fvc::flux(U)
+);
+
+if (args.found("initialiseUBCs"))
+{
+    U.correctBoundaryConditions();
+    phi = fvc::flux(U);
+}
+
+
+// Construct a pressure field
+// If it is available read it otherwise construct from the velocity BCs
+// converting fixed-value BCs to zero-gradient and vice versa.
+
+// Allow override from command-line -pName option
+const word pName = args.getOrDefault<word>("pName", "p");
+
+// Infer the pressure BCs from the velocity
+wordList pBCTypes
+(
+    U.boundaryField().size(),
+    fixedValueFvPatchScalarField::typeName
+);
+
+forAll(U.boundaryField(), patchi)
+{
+    if (U.boundaryField()[patchi].fixesValue())
+    {
+        pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
+    }
+}
+
+// Note that registerObject is false for the pressure field. The pressure
+// field in this solver doesn't have a physical value during the solution.
+// It shouldn't be looked up and used by sub models or boundary conditions.
+Info<< "Constructing pressure field " << pName << nl << endl;
+volScalarField p
+(
+    IOobject
+    (
+        pName,
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE,
+        false
+    ),
+    mesh,
+    dimensionedScalar(sqr(dimVelocity), Zero),
+    pBCTypes
+);
+
+// Infer the velocity potential BCs from the pressure
+wordList PhiBCTypes
+(
+    p.boundaryField().size(),
+    zeroGradientFvPatchScalarField::typeName
+);
+
+forAll(p.boundaryField(), patchi)
+{
+    if (p.boundaryField()[patchi].fixesValue())
+    {
+        PhiBCTypes[patchi] = fixedValueFvPatchScalarField::typeName;
+    }
+}
+
+Info<< "Constructing velocity potential field Phi\n" << endl;
+volScalarField Phi
+(
+    IOobject
+    (
+        "Phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimLength*dimVelocity, Zero),
+    PhiBCTypes
+);
+
+label PhiRefCell = 0;
+scalar PhiRefValue = 0;
+setRefCell
+(
+    Phi,
+    potentialFlow.dict(),
+    PhiRefCell,
+    PhiRefValue
+);
+mesh.setFluxRequired(Phi.name());
+
+#include "createMRF.H"

 // Add solver-specific interpolations
 {
--- 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.
@ -83,7 +83,6 @@ Description
    \heading Options
    \plaintable
        -writep   | write the Euler pressure
-        -writephi | Write the final volumetric flux
        -writePhi | Write the final velocity potential
        -initialiseUBCs | Update the velocity boundaries before solving for Phi
    \endplaintable
@ -118,12 +117,6 @@ int main(int argc, char *argv[])
        "Initialise U boundary conditions"
    );

-    argList::addBoolOption
-    (
-        "writephi",
-        "Write the final volumetric flux field"
-    );
-
    argList::addBoolOption
    (
        "writePhi",
@ -142,8 +135,6 @@ int main(int argc, char *argv[])
        "Execute functionObjects"
    );

-    #include "addRegionOption.H"
-    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createNamedDynamicFvMesh.H"
@ -158,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
@ -203,16 +195,11 @@ int main(int argc, char *argv[])
            << endl;
    }

-    // Write U
+    // Write U and phi
    U.write();
+    phi.write();

-    // Optionally write the volumetric flux, phi
-    if (args.found("writephi"))
-    {
-        phi.write();
-    }
-
-    // Optionally write velocity potential, Phi
+    // Optionally write Phi
    if (args.found("writePhi"))
    {
        Phi.write();
--- 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/PDRFoam.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoam.C
@ -160,7 +160,7 @@ int main(int argc, char *argv[])
        runTime.printExecutionTime(Info);
    }

-    Info<< "End\n";
+    Info<< "\n end\n";

    return 0;
 }
--- 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;
@ -207,7 +207,7 @@ int main(int argc, char *argv[])
        runTime.printExecutionTime(Info);
    }

-    Info<< "End\n";
+    Info<< "\n end\n";

    return 0;
 }
--- a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -43,9 +43,9 @@ Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New

    Info<< "Selecting drag model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -57,9 +57,7 @@ Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New
    }

    return autoPtr<PDRDragModel>
-    (
-        ctorPtr(dict, turbulence, rho, U, phi)
-    );
+        (cstrIter()(dict, turbulence, rho, U, phi));
 }


--- a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -94,14 +93,24 @@ Foam::PDRDragModels::basic::~basic()

 Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basic::Dcu() const
 {
-    auto tDragDcu = volSymmTensorField::New
+    tmp<volSymmTensorField> tDragDcu
    (
-        "tDragDcu",
-        IOobject::NO_REGISTER,
-        U_.mesh(),
-        dimensionedSymmTensor(dimMass/dimTime/dimVolume, Zero)
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                "tDragDcu",
+                U_.mesh().time().constant(),
+                U_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U_.mesh(),
+            dimensionedSymmTensor(dimMass/dimTime/dimVolume, Zero)
+        )
    );
-    auto& DragDcu = tDragDcu.ref();
+
+    volSymmTensorField& DragDcu = tDragDcu.ref();

    if (on_)
    {
@ -118,14 +127,24 @@ Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basic::Dcu() const

 Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basic::Gk() const
 {
-    auto tGk = volScalarField::New
+    tmp<volScalarField> tGk
    (
-        "tGk",
-        IOobject::NO_REGISTER,
-        U_.mesh(),
-        dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
+        new volScalarField
+        (
+            IOobject
+            (
+                "tGk",
+                U_.mesh().time().constant(),
+                U_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U_.mesh(),
+            dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
+        )
    );
-    auto& Gk = tGk.ref();
+
+    volScalarField& Gk = tGk.ref();

    if (on_)
    {
--- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020,2023 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -132,11 +132,6 @@ void PDRkEpsilon::correct()

    // Update epsilon and G at the wall
    epsilon_.boundaryFieldRef().updateCoeffs();
-    // Push new cell values to
-    // coupled neighbours. Note that we want to avoid the re-updateCoeffs
-    // of the wallFunctions so make sure to bypass the evaluate on
-    // those patches and only do the coupled ones.
-    epsilon_.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();

    // Add the blockage generation term so that it is included consistently
    // in both the k and epsilon equations
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C
@ -69,8 +69,7 @@ Foam::XiEqModels::Gulder::~Gulder()
 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 {
    volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
-    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
-    const volScalarField& epsilon = tepsilon();
+    const volScalarField& epsilon = turbulence_.epsilon();

    if (subGridSchelkin_)
    {
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -76,10 +75,8 @@ Foam::XiEqModels::SCOPEXiEq::~SCOPEXiEq()

 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
 {
-    const tmp<volScalarField> tk(turbulence_.k());
-    const volScalarField& k = tk();
-    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
-    const volScalarField& epsilon = tepsilon();
+    const volScalarField& k = turbulence_.k();
+    const volScalarField& epsilon = turbulence_.epsilon();

    volScalarField up(sqrt((2.0/3.0)*k));
    if (subGridSchelkin_)
@ -94,14 +91,23 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
    volScalarField K(0.157*upBySu/sqrt(Rl));
    volScalarField Ma(MaModel.Ma());

-    auto tXiEq = volScalarField::New
+    tmp<volScalarField> tXiEq
    (
-        "XiEq",
-        IOobject::NO_REGISTER,
-        epsilon.mesh(),
-        dimensionedScalar(dimless, Zero)
+        new volScalarField
+        (
+            IOobject
+            (
+                "XiEq",
+                epsilon.time().timeName(),
+                epsilon.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            epsilon.mesh(),
+            dimensionedScalar(dimless, Zero)
+        )
    );
-    auto& xieq = tXiEq.ref();
+    volScalarField& xieq = tXiEq.ref();

    forAll(xieq, celli)
    {
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C
@ -78,10 +78,11 @@ bool Foam::XiEqModel::read(const dictionary& XiEqProperties)
 void Foam::XiEqModel::writeFields() const
 {
    //***HGW It is not clear why B is written here
-    const auto* B = Su_.mesh().cfindObject<volSymmTensorField>("B");
-    if (B)
+    if (Su_.mesh().foundObject<volSymmTensorField>("B"))
    {
-        B->write();
+        const volSymmTensorField& B =
+            Su_.mesh().lookupObject<volSymmTensorField>("B");
+        B.write();
    }
 }

@ -97,26 +98,39 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
    const volSymmTensorField& nsv =
        mesh.lookupObject<volSymmTensorField>("nsv");

-    auto tN = volScalarField::New
+    tmp<volScalarField> tN
    (
-        "tN",
-        IOobject::NO_REGISTER,
-        mesh,
-        dimensionedScalar(Nv.dimensions(), Zero)
+        new volScalarField
+        (
+            IOobject
+            (
+                "tN",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar(Nv.dimensions(), Zero)
+        )
    );
-    auto& N = tN.ref();
-
+    volScalarField& N = tN.ref();
    N.primitiveFieldRef() = Nv.primitiveField()*pow(mesh.V(), 2.0/3.0);

-    auto tns = volSymmTensorField::New
+    volSymmTensorField ns
    (
-        "tns",
-        IOobject::NO_REGISTER,
+        IOobject
+        (
+            "tns",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
        mesh,
        dimensionedSymmTensor(nsv.dimensions(), Zero)
    );
-    auto& ns = tns.ref();
-
    ns.primitiveFieldRef() = nsv.primitiveField()*pow(mesh.V(), 2.0/3.0);

    const volVectorField Uhat
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModelNew.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -42,9 +42,9 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New

    Info<< "Selecting flame-wrinkling model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -55,7 +55,7 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
        ) << exit(FatalIOError);
    }

-    return autoPtr<XiEqModel>(ctorPtr(dict, thermo, turbulence, Su));
+    return autoPtr<XiEqModel>(cstrIter()(dict, thermo, turbulence, Su));
 }


--- a/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C
@ -66,8 +66,7 @@ Foam::XiGModels::KTS::~KTS()
 Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::G() const
 {
    volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
-    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
-    const volScalarField& epsilon = tepsilon();
+    const volScalarField& epsilon = turbulence_.epsilon();

    volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));

--- a/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModelNew.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -42,9 +42,9 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New

    Info<< "Selecting flame-wrinkling model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -55,7 +55,7 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
        ) << exit(FatalIOError);
    }

-    return autoPtr<XiGModel>(ctorPtr(dict, thermo, turbulence, Su));
+    return autoPtr<XiGModel>(cstrIter()(dict, thermo, turbulence, Su));
 }


--- a/applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModelNew.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModelNew.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -45,9 +45,9 @@ Foam::autoPtr<Foam::XiModel> Foam::XiModel::New

    Info<< "Selecting flame-wrinkling model " << modelType << endl;

-    auto* ctorPtr = dictionaryConstructorTable(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

-    if (!ctorPtr)
+    if (!cstrIter.found())
    {
        FatalIOErrorInLookup
        (
@ -59,9 +59,7 @@ Foam::autoPtr<Foam::XiModel> Foam::XiModel::New
    }

    return autoPtr<XiModel>
-    (
-        ctorPtr(dict, thermo, turbulence, Su, rho, b, phi)
-    );
+        (cstrIter()(dict, thermo, turbulence, Su, rho, b, phi));
 }


--- a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
+++ b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -256,14 +255,24 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
    scalar phi
 ) const
 {
-    auto tSu0 = volScalarField::New
+    tmp<volScalarField> tSu0
    (
-        "Su0",
-        IOobject::NO_REGISTER,
-        p.mesh(),
-        dimensionedScalar(dimVelocity, Zero)
+        new volScalarField
+        (
+            IOobject
+            (
+                "Su0",
+                p.time().timeName(),
+                p.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            p.mesh(),
+            dimensionedScalar(dimVelocity, Zero)
+        )
    );
-    auto& Su0 = tSu0.ref();
+
+    volScalarField& Su0 = tSu0.ref();

    forAll(Su0, celli)
    {
@ -295,14 +304,24 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
    const volScalarField& phi
 ) const
 {
-    auto tSu0 = volScalarField::New
+    tmp<volScalarField> tSu0
    (
-        "Su0",
-        IOobject::NO_REGISTER,
-        p.mesh(),
-        dimensionedScalar(dimVelocity, Zero)
+        new volScalarField
+        (
+            IOobject
+            (
+                "Su0",
+                p.time().timeName(),
+                p.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            p.mesh(),
+            dimensionedScalar(dimVelocity, Zero)
+        )
    );
-    auto& Su0 = tSu0.ref();
+
+    volScalarField& Su0 = tSu0.ref();

    forAll(Su0, celli)
    {
@ -339,14 +358,24 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
    const volScalarField& phi
 ) const
 {
-    auto tMa = volScalarField::New
+    tmp<volScalarField> tMa
    (
-        "Ma",
-        IOobject::NO_REGISTER,
-        phi.mesh(),
-        dimensionedScalar(dimless, Zero)
+        new volScalarField
+        (
+            IOobject
+            (
+                "Ma",
+                phi.time().timeName(),
+                phi.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            phi.mesh(),
+            dimensionedScalar(dimless, Zero)
+        )
    );
-    auto& ma = tMa.ref();
+
+    volScalarField& ma = tMa.ref();

    forAll(ma, celli)
    {
@ -389,12 +418,21 @@ Foam::laminarFlameSpeedModels::SCOPE::Ma() const
    {
        const fvMesh& mesh = psiuReactionThermo_.p().mesh();

-        return volScalarField::New
+        return tmp<volScalarField>
        (
-            "Ma",
-            IOobject::NO_REGISTER,
-            mesh,
-            dimensionedScalar("Ma", dimless, Ma(equivalenceRatio_))
+            new volScalarField
+            (
+                IOobject
+                (
+                    "Ma",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("Ma", dimless, Ma(equivalenceRatio_))
+            )
        );
    }
 }
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@ -27,7 +27,7 @@ if (pimple.transonic())
            betav*fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -57,7 +57,7 @@ else
            betav*fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
--- 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/XiDyMFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H
@ -30,7 +30,7 @@ if (pimple.transonic())
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -66,7 +66,7 @@ else
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
--- 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/XiFoam/XiEngineFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/XiEngineFoam/pEqn.H
@ -35,7 +35,7 @@ if (pimple.transonic())
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -71,7 +71,7 @@ else
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
--- a/applications/solvers/combustion/XiFoam/bEqn.H
+++ b/applications/solvers/combustion/XiFoam/bEqn.H
@ -171,7 +171,10 @@ if (ign.ignited())

        fvOptions.correct(Su);

-        Su.clamp_range(SuMin, SuMax);
+        // Limit the maximum Su
+        // ~~~~~~~~~~~~~~~~~~~~
+        Su.min(SuMax);
+        Su.max(SuMin);
    }
    else
    {
@ -215,7 +218,7 @@ if (ign.ignited())
          + (
                scalar(1)
              + (2*XiShapeCoef)
-               *(scalar(0.5) - clamp(b, zero_one{}))
+               *(scalar(0.5) - min(max(b, scalar(0)), scalar(1)))
            )*(XiEqStar - scalar(1.001))
        );

@ -223,7 +226,7 @@ if (ign.ignited())
        volScalarField R(Gstar*XiEqStar/(XiEqStar - scalar(1)));
        volScalarField G(R*(XiEq - scalar(1.001))/XiEq);

-        //R *= (Gstar + 2*mag(devSymm(fvc::grad(U))))/Gstar;
+        //R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;

        // Solve for the flame wrinkling
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@ -29,7 +29,7 @@ if (pimple.transonic())
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -64,7 +64,7 @@ else
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
--- a/applications/solvers/combustion/chemFoam/Make/options
+++ b/applications/solvers/combustion/chemFoam/Make/options
@ -6,13 +6,11 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
-    -lchemistryModel \
-    -lthermophysicalProperties
+    -lchemistryModel
--- a/applications/solvers/combustion/chemFoam/YEqn.H
+++ b/applications/solvers/combustion/chemFoam/YEqn.H
@ -6,7 +6,7 @@
        solve
        (
            fvm::ddt(rho, Yi) - chemistry.RR(specieI),
-            "Yi"
+            mesh.solver("Yi")
        );
    }
 }
--- a/applications/solvers/combustion/chemFoam/createBaseFields.H
+++ b/applications/solvers/combustion/chemFoam/createBaseFields.H
@ -12,7 +12,7 @@ Info<< "Creating base fields for time " << runTime.timeName() << endl;
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
-            IOobject::NO_REGISTER
+            false
        ),
        mesh,
        dimensionedScalar("Ydefault", dimless, 1)
@ -29,7 +29,7 @@ Info<< "Creating base fields for time " << runTime.timeName() << endl;
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
-            IOobject::NO_REGISTER
+            false
        ),
        mesh,
        dimensionedScalar("p", dimPressure, p0)
@ -46,7 +46,7 @@ Info<< "Creating base fields for time " << runTime.timeName() << endl;
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
-            IOobject::NO_REGISTER
+            false
        ),
        mesh,
        dimensionedScalar("T", dimTemperature, T0)
--- a/applications/solvers/combustion/chemFoam/setDeltaT.H
+++ b/applications/solvers/combustion/chemFoam/setDeltaT.H
@ -1,5 +1,5 @@
 if (adjustTimeStep)
 {
    runTime.setDeltaT(min(dtChem, maxDeltaT));
-    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+    Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
 }
--- a/applications/solvers/combustion/chemFoam/solveChemistry.H
+++ b/applications/solvers/combustion/chemFoam/solveChemistry.H
@ -1,3 +1,3 @@
-dtChem = chemistry.solve(runTime.deltaTValue());
+dtChem = chemistry.solve(runTime.deltaT().value());
 scalar Qdot = chemistry.Qdot()()[0]/rho[0];
-integratedHeat += Qdot*runTime.deltaTValue();
+integratedHeat += Qdot*runTime.deltaT().value();
--- 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,7 @@
 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 \
@ -15,6 +14,7 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
@ -23,9 +23,7 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/pyrolysisModels/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/faOptions/lnInclude
+    -I$(LIB_SRC)/ODE/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -34,7 +32,6 @@ EXE_LIBS = \
    -lsampling \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
@ -53,7 +50,4 @@ EXE_LIBS = \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -lODE \
-    -lregionFaModels \
-    -lfiniteArea \
-    -lfaOptions
+    -lODE
--- a/applications/solvers/combustion/fireFoam/YEEqn.H
+++ b/applications/solvers/combustion/fireFoam/YEEqn.H
@ -35,17 +35,17 @@ tmp<fv::convectionScheme<scalar>> mvConvection

            fvOptions.constrain(YiEqn);

-            YiEqn.solve("Yi");
+            YiEqn.solve(mesh.solver("Yi"));

            fvOptions.correct(Yi);

-            Yi.clamp_min(0);
+            Yi.max(0.0);
            Yt += Yi;
        }
    }

    Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].clamp_min(0);
+    Y[inertIndex].max(0.0);

    radiation->correct();

--- 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/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@ -36,7 +36,7 @@ while (pimple.correctNonOrthogonal())
      + fvOptions(psi, p_rgh, rho.name())
    );

-    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

    if (pimple.finalNonOrthogonalIter())
    {
--- a/applications/solvers/combustion/fireFoam/phrghEqn.H
+++ b/applications/solvers/combustion/fireFoam/phrghEqn.H
@ -10,8 +10,7 @@ if (pimple.dict().getOrDefault("hydrostaticInitialization", false))
                "0",
                mesh,
                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                IOobject::REGISTER
+                IOobject::NO_WRITE
            ),
            mesh
        )
--- 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/YEqn.H
+++ b/applications/solvers/combustion/reactingFoam/YEqn.H
@ -34,15 +34,15 @@ tmp<fv::convectionScheme<scalar>> mvConvection

            fvOptions.constrain(YiEqn);

-            YiEqn.solve("Yi");
+            YiEqn.solve(mesh.solver("Yi"));

            fvOptions.correct(Yi);

-            Yi.clamp_min(0);
+            Yi.max(0.0);
            Yt += Yi;
        }
    }

    Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].clamp_min(0);
+    Y[inertIndex].max(0.0);
 }
--- 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/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@ -37,7 +37,7 @@ if (pimple.transonic())
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -72,7 +72,7 @@ else
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
--- a/applications/solvers/combustion/reactingFoam/pcEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pcEqn.H
@ -49,7 +49,7 @@ if (pimple.transonic())
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -89,7 +89,7 @@ else
            fvOptions(psi, p, rho.name())
        );

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
--- 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/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@ -41,7 +41,7 @@ while (pimple.correctNonOrthogonal())
      - fvm::laplacian(rhorAUf, p_rgh)
    );

-    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

    if (pimple.finalNonOrthogonalIter())
    {
--- 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/BCs/T/smoluchowskiJumpTFvPatchScalarField.C
+++ b/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C
@ -103,10 +103,16 @@ Foam::smoluchowskiJumpTFvPatchScalarField::smoluchowskiJumpTFvPatchScalarField
            << exit(FatalIOError);
    }

-    if (!this->readValueEntry(dict))
+    if (dict.found("value"))
    {
-        // Fallback: set to the internal field
-        fvPatchField<scalar>::patchInternalField(*this);
+        fvPatchField<scalar>::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchField<scalar>::operator=(patchInternalField());
    }

    refValue() = *this;
@ -159,10 +165,14 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::updateCoeffs()
        return;
    }

-    const auto& pmu = patch().lookupPatchField<volScalarField>(muName_);
-    const auto& prho = patch().lookupPatchField<volScalarField>(rhoName_);
-    const auto& ppsi = patch().lookupPatchField<volScalarField>(psiName_);
-    const auto& pU = patch().lookupPatchField<volVectorField>(UName_);
+    const fvPatchScalarField& pmu =
+        patch().lookupPatchField<volScalarField, scalar>(muName_);
+    const fvPatchScalarField& prho =
+        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+    const fvPatchField<scalar>& ppsi =
+        patch().lookupPatchField<volScalarField, scalar>(psiName_);
+    const fvPatchVectorField& pU =
+        patch().lookupPatchField<volVectorField, vector>(UName_);

    // Prandtl number reading consistent with rhoCentralFoam
    const dictionary& thermophysicalProperties =
@ -197,7 +207,7 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::updateCoeffs()
 // Write
 void Foam::smoluchowskiJumpTFvPatchScalarField::write(Ostream& os) const
 {
-    fvPatchField<scalar>::write(os);
+    fvPatchScalarField::write(os);

    os.writeEntryIfDifferent<word>("U", "U", UName_);
    os.writeEntryIfDifferent<word>("rho", "rho", rhoName_);
@ -207,7 +217,7 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::write(Ostream& os) const
    os.writeEntry("accommodationCoeff", accommodationCoeff_);
    Twall_.writeEntry("Twall", os);
    os.writeEntry("gamma", gamma_);
-    fvPatchField<scalar>::writeValueEntry(os);
+    writeEntry("value", os);
 }


--- a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.H
+++ b/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.H
@ -110,6 +110,15 @@ public:
            const fvPatchFieldMapper&
        );

+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new smoluchowskiJumpTFvPatchScalarField(*this)
+            );
+        }
+
        //- Construct as copy setting internal field reference
        smoluchowskiJumpTFvPatchScalarField
        (
@ -117,22 +126,18 @@ public:
            const DimensionedField<scalar, volMesh>&
        );

-        //- Return a clone
-        virtual tmp<fvPatchField<scalar>> clone() const
-        {
-            return fvPatchField<scalar>::Clone(*this);
-        }
-
-        //- Clone with an internal field reference
-        virtual tmp<fvPatchField<scalar>> clone
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
        (
            const DimensionedField<scalar, volMesh>& iF
        ) const
        {
-            return fvPatchField<scalar>::Clone(*this, iF);
+            return tmp<fvPatchScalarField>
+            (
+                new smoluchowskiJumpTFvPatchScalarField(*this, iF)
+            );
        }

-
        // Mapping functions

            //- Map (and resize as needed) from self given a mapping object
--- a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
+++ b/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
@ -105,15 +105,18 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
            << exit(FatalIOError);
    }

-    if (this->readValueEntry(dict))
+    if (dict.found("value"))
    {
-        const auto* hasRefValue = dict.findEntry("refValue", keyType::LITERAL);
-        const auto* hasFrac = dict.findEntry("valueFraction", keyType::LITERAL);
+        fvPatchField<vector>::operator=
+        (
+            vectorField("value", dict, p.size())
+        );

-        if (hasRefValue && hasFrac)
+        if (dict.found("refValue") && dict.found("valueFraction"))
        {
-            this->refValue().assign(*hasRefValue, p.size());
-            this->valueFraction().assign(*hasFrac, p.size());
+            this->refValue() = vectorField("refValue", dict, p.size());
+            this->valueFraction() =
+                scalarField("valueFraction", dict, p.size());
        }
        else
        {
@ -152,9 +155,12 @@ void Foam::maxwellSlipUFvPatchVectorField::updateCoeffs()
        return;
    }

-    const auto& pmu = patch().lookupPatchField<volScalarField>(muName_);
-    const auto& prho = patch().lookupPatchField<volScalarField>(rhoName_);
-    const auto& ppsi = patch().lookupPatchField<volScalarField>(psiName_);
+    const fvPatchScalarField& pmu =
+        patch().lookupPatchField<volScalarField, scalar>(muName_);
+    const fvPatchScalarField& prho =
+        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+    const fvPatchField<scalar>& ppsi =
+        patch().lookupPatchField<volScalarField, scalar>(psiName_);

    Field<scalar> C1
    (
@ -181,8 +187,8 @@ void Foam::maxwellSlipUFvPatchVectorField::updateCoeffs()

    if (curvature_)
    {
-        const auto& ptauMC =
-            patch().lookupPatchField<volTensorField>(tauMCName_);
+        const fvPatchTensorField& ptauMC =
+            patch().lookupPatchField<volTensorField, tensor>(tauMCName_);
        vectorField n(patch().nf());

        refValue() -= C1/prho*transform(I - n*n, (n & ptauMC));
@ -194,7 +200,7 @@ void Foam::maxwellSlipUFvPatchVectorField::updateCoeffs()

 void Foam::maxwellSlipUFvPatchVectorField::write(Ostream& os) const
 {
-    fvPatchField<vector>::write(os);
+    fvPatchVectorField::write(os);
    os.writeEntryIfDifferent<word>("T", "T", TName_);
    os.writeEntryIfDifferent<word>("rho", "rho", rhoName_);
    os.writeEntryIfDifferent<word>("psi", "thermo:psi", psiName_);
@ -209,7 +215,7 @@ void Foam::maxwellSlipUFvPatchVectorField::write(Ostream& os) const
    refValue().writeEntry("refValue", os);
    valueFraction().writeEntry("valueFraction", os);

-    fvPatchField<vector>::writeValueEntry(os);
+    writeEntry("value", os);
 }


--- a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.H
+++ b/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.H
@ -118,6 +118,15 @@ public:
            const fvPatchFieldMapper&
        );

+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new maxwellSlipUFvPatchVectorField(*this)
+            );
+        }
+
        //- Construct as copy setting internal field reference
        maxwellSlipUFvPatchVectorField
        (
@ -125,19 +134,16 @@ public:
            const DimensionedField<vector, volMesh>&
        );

-        //- Return a clone
-        virtual tmp<fvPatchField<vector>> clone() const
-        {
-            return fvPatchField<vector>::Clone(*this);
-        }
-
-        //- Clone with an internal field reference
-        virtual tmp<fvPatchField<vector>> clone
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
        (
            const DimensionedField<vector, volMesh>& iF
        ) const
        {
-            return fvPatchField<vector>::Clone(*this, iF);
+            return tmp<fvPatchVectorField>
+            (
+                new maxwellSlipUFvPatchVectorField(*this, iF)
+            );
        }


--- a/applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.C
+++ b/applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.C
@ -104,8 +104,11 @@ void Foam::fixedRhoFvPatchScalarField::updateCoeffs()
        return;
    }

-    const auto& psip = patch().lookupPatchField<volScalarField>(psiName_);
-    const auto& pp = patch().lookupPatchField<volScalarField>(pName_);
+    const fvPatchField<scalar>& psip =
+        patch().lookupPatchField<volScalarField, scalar>(psiName_);
+
+    const fvPatchField<scalar>& pp =
+        patch().lookupPatchField<volScalarField, scalar>(pName_);

    operator==(psip*pp);

@ -115,10 +118,11 @@ void Foam::fixedRhoFvPatchScalarField::updateCoeffs()

 void Foam::fixedRhoFvPatchScalarField::write(Ostream& os) const
 {
-    fvPatchField<scalar>::write(os);
+    fvPatchScalarField::write(os);
+
    os.writeEntryIfDifferent<word>("p", "p", pName_);
    os.writeEntryIfDifferent<word>("psi", "thermo:psi", psiName_);
-    fvPatchField<scalar>::writeValueEntry(os);
+    writeEntry("value", os);
 }

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
--- a/applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.H
+++ b/applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.H
@ -135,6 +135,15 @@ public:
            const fixedRhoFvPatchScalarField&
        );

+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new fixedRhoFvPatchScalarField(*this)
+            );
+        }
+
        //- Construct as copy setting internal field reference
        fixedRhoFvPatchScalarField
        (
@ -142,19 +151,16 @@ public:
            const DimensionedField<scalar, volMesh>&
        );

-        //- Return a clone
-        virtual tmp<fvPatchField<scalar>> clone() const
-        {
-            return fvPatchField<scalar>::Clone(*this);
-        }
-
-        //- Clone with an internal field reference
-        virtual tmp<fvPatchField<scalar>> clone
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
        (
            const DimensionedField<scalar, volMesh>& iF
        ) const
        {
-            return fvPatchField<scalar>::Clone(*this, iF);
+            return tmp<fvPatchScalarField>
+            (
+                new fixedRhoFvPatchScalarField(*this, iF)
+            );
        }


--- 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/createFieldRefs.H
+++ b/applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H
@ -1,9 +1,10 @@
 volScalarField& p = thermo.p();
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
+const volScalarField& mu = thermo.mu();

 bool inviscid(true);
-if (max(thermo.mu().cref().primitiveField()) > 0.0)
+if (max(mu.primitiveField()) > 0.0)
 {
    inviscid = false;
 }
--- a/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H
+++ b/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H
@ -17,7 +17,7 @@ tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
            vf,
            dir,
            "reconstruct("
-          + (reconFieldName.empty() ? vf.name() : reconFieldName)
+          + (reconFieldName != word::null ? reconFieldName : vf.name())
          + ')'
        )
    );
--- 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/rhoPimpleAdiabaticFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleAdiabaticFoam/pEqn.H
@ -63,7 +63,7 @@
                fvOptions(psi, p, rho.name())
            );

-            pEqn.solve(p.select(pimple.finalInnerIter()));
+            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

            // Rhie & Chow interpolation (part 2)
            if (pimple.finalNonOrthogonalIter())
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@ -8,8 +8,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)/regionFaModels/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -20,11 +19,8 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -lsampling \
-    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -latmosphericModels
--- 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/correctPhi.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/correctPhi.H
@ -29,7 +29,7 @@ if (mesh.changing())
    wordList pcorrTypes
    (
        p.boundaryField().size(),
-        fvPatchFieldBase::zeroGradientType()
+        zeroGradientFvPatchScalarField::typeName
    );

    // Set BCs of pcorr to fixed-value for patches at which p is fixed
@ -72,7 +72,7 @@ if (mesh.changing())
            divrhoU()
        );

-        pcorrEqn.solve(pcorr.select(pimple.finalInnerIter()));
+        pcorrEqn.solve(mesh.solver(pcorr.select(pimple.finalInnerIter())));
        //Bypass virtual layer
        //mesh.fvMesh::solve(pcorrEqn, d);

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

@ -56,7 +67,7 @@ if (pimple.transonic())
        // Relax the pressure equation to ensure diagonal-dominance
        pEqn.relax();

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -78,7 +89,7 @@ else
    {
        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -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/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@ -54,7 +54,7 @@ if (pimple.transonic())
        // Relax the pressure equation to ensure diagonal-dominance
        pEqn.relax();

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -76,7 +76,7 @@ else
    {
        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@ -65,7 +65,7 @@ if (pimple.transonic())
        // Relax the pressure equation to ensure diagonal-dominance
        pEqn.relax();

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
@ -90,7 +90,7 @@ else
    {
        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));

-        pEqn.solve(p.select(pimple.finalInnerIter()));
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
--- 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
@ -7,7 +7,6 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
@ -19,7 +18,4 @@ EXE_LIBS = \
    -lspecie \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
-    -lthermoTools \
-    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -latmosphericModels
--- a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/EEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/EEqn.H
@ -0,0 +1,26 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+}
--- a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/files
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/files
@ -0,0 +1,3 @@
+SRFrhoSimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/SRFrhoSimpleFoam
--- a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/options
@ -0,0 +1,21 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -latmosphericModels
--- a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/SRFrhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/SRFrhoSimpleFoam.C
@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    SRFSimpleFoam
+
+Group
+    grpIncompressibleSolvers
+
+Description
+    Steady-state solver for incompressible, turbulent flow of non-Newtonian
+    fluids in a single rotating frame.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "fluidThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "simpleControl.H"
+#include "pressureControl.H"
+#include "fvOptions.H"
+#include "SRFModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Steady-state solver for compressible, turbulent flow"
+        " of non-Newtonian fluids in a single rotating frame."
+    );
+
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "initContinuityErrs.H"
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity SIMPLE corrector
+        {
+            #include "UrelEqn.H"
+            #include "EEqn.H"
+            #include "pEqn.H"
+        }
+
+        U = Urel + SRF->U();
+
+        turbulence->correct();
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/UrelEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/UrelEqn.H
@ -0,0 +1,22 @@
+    // Relative momentum predictor
+
+    tmp<fvVectorMatrix> tUrelEqn
+    (
+        fvm::div(phi, Urel)
+      + turbulence->divDevRhoReff(Urel)
+      + rho*SRF->Su()
+     ==
+        fvOptions(rho, Urel)
+    );
+    fvVectorMatrix& UrelEqn = tUrelEqn.ref();
+
+    UrelEqn.relax();
+
+    fvOptions.constrain(UrelEqn);
+
+    if (simple.momentumPredictor())
+    {
+        solve(UrelEqn == -fvc::grad(p));
+
+        fvOptions.correct(Urel);
+    }
--- a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFieldRefs.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFieldRefs.H
@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();
--- a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFields.H
@ -0,0 +1,89 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+autoPtr<fluidThermo> pThermo
+(
+    fluidThermo::New(mesh)
+);
+fluidThermo& thermo = pThermo();
+thermo.validate(args.executable(), "h", "e");
+
+volScalarField& p = thermo.p();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field Urel\n" << endl;
+volVectorField Urel
+(
+    IOobject
+    (
+        "Urel",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(rho*Urel) & mesh.Sf()
+);
+
+pressureControl pressureControl(p, rho, simple.dict());
+
+mesh.setFluxRequired(p.name());
+
+Info<< "Creating SRF model\n" << endl;
+autoPtr<SRF::SRFModel> SRF(SRF::SRFModel::New(Urel));
+
+// Construct the absolute velocity
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    Urel + SRF->U()
+);
+
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+#include "createFvOptions.H"
--- a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/pEqn.H
@ -0,0 +1,109 @@
+volScalarField rAUrel(1.0/UrelEqn.A());
+surfaceScalarField rhorAUrelf("rhorAUf", fvc::interpolate(rho*rAUrel));
+volVectorField HbyA(constrainHbyA(rAUrel*UrelEqn.H(), Urel, p));
+tUrelEqn.clear();
+
+bool closedVolume = false;
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, Urel, phiHbyA, rhorAUrelf);
+
+if (simple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+    );
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rhorAUrelf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+else
+{
+    closedVolume = adjustPhi(phiHbyA, Urel, p);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUrelf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "incompressible/continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+Urel = HbyA - rAUrel*fvc::grad(p);
+Urel.correctBoundaryConditions();
+fvOptions.correct(Urel);
+
+bool pLimited = pressureControl.limit(p);
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+}
+
+if (pLimited || closedVolume)
+{
+    p.correctBoundaryConditions();
+}
+
+rho = thermo.rho();
+
+if (!simple.transonic())
+{
+    rho.relax();
+}
--- 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
+            (
+                dynamicFvMesh::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/Show More
+++ b/Show More
Author	SHA1	Message	Date
Andrew Heather	a62529c0e9	ENH: new SRFrhoSimpleFoam solver	2021-03-16 21:36:43 +00:00
Andrew Heather	9423976dda	ENH: SRF - refactoring and enabled ramping of rotation	2021-03-16 21:36:43 +00:00