TUT: added two unsteady adjoint optimisation tutorials

showcasing the use of compressedFullStorage and binomialCheckPointing

.gitlab/issue_templates/bug.md

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

.gitmodules

@@ -1,8 +1,8 @@
 [submodule "cfmesh"]
-	path = plugins/cfmesh
+	path = modules/cfmesh
 	url = https://develop.openfoam.com/Community/integration-cfmesh.git
 [submodule "avalanche"]
-	path = plugins/avalanche
+	path = modules/avalanche
 	url = https://develop.openfoam.com/Community/avalanche.git
 [submodule "adios"]
 	path = modules/adios
@@ -17,9 +17,3 @@
 [submodule "external-solver"]
 	path = modules/external-solver
 	url = https://develop.openfoam.com/Modules/external-solver.git
-[submodule "turbulence-community"]
-	path = plugins/turbulence-community
-	url = https://gitlab.com/openfoam/community/tc-turbulence/turbulence-community.git
-[submodule "plugins/data-community"]
-	path = plugins/data-community
-	url = https://gitlab.com/openfoam/community/sig-data-modelling/data-community.git

Allwmake

@@ -68,34 +68,19 @@ src/Allwmake $targetType $*
 # OpenFOAM applications
 applications/Allwmake $targetType $*
 
-#------------------------------------------------------------------------------
-# Additional components
 
+# Additional components/modules
 case "$FOAM_MODULE_PREFIX" in
 (false | none)
     echo ========================================
     echo "OpenFOAM modules disabled (prefix=${FOAM_MODULE_PREFIX})"
-    echo "Can be built separately:"
-    echo
-    echo "    ./Allwmake-modules -prefix=..."
-    echo
-    echo ========================================
     echo
     ;;
 (*)
     # Use wmake -all instead of Allwmake to allow for overrides
     ( cd "$WM_PROJECT_DIR/modules" 2>/dev/null && wmake -all )
-
-    echo ========================================
-    echo "The optional plugins can be built separately:"
-    echo
-    echo "    ./Allwmake-plugins -prefix=..."
-    echo
-    echo ========================================
-    echo
 esac
 
-#------------------------------------------------------------------------------
 # Count files in given directory. Ignore "Test-*" binaries.
 _foamCountDirEntries()
 {
@@ -111,12 +96,10 @@ 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 ========================================
 

Allwmake-modules

@@ -1,37 +0,0 @@
-#!/bin/sh
-# Run from OPENFOAM top-level directory only
-cd "${0%/*}" || exit
-wmake -check-dir "$WM_PROJECT_DIR" 2>/dev/null || {
-    echo "Error (${0##*/}) : not located in \$WM_PROJECT_DIR"
-    echo "    Check your OpenFOAM environment and installation"
-    exit 1
-}
-if [ -f "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments ]
-then  . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments -no-recursion "$@" || \
-    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
-
-#------------------------------------------------------------------------------
-# Additional components
-
-case "$FOAM_MODULE_PREFIX" in
-(false | none)
-    echo ========================================
-    echo "OpenFOAM modules disabled (prefix=${FOAM_MODULE_PREFIX})"
-    echo "Can be built separately:"
-    echo
-    echo "    ./Allwmake-modules -prefix=..."
-    echo
-    echo ========================================
-    echo
-    ;;
-(*)
-    # Use wmake -all instead of Allwmake to allow for overrides
-    ( cd "$WM_PROJECT_DIR/modules" 2>/dev/null && wmake -all )
-esac
-
-#------------------------------------------------------------------------------

Allwmake-plugins

@@ -1,37 +0,0 @@
-#!/bin/sh
-# Run from OPENFOAM top-level directory only
-cd "${0%/*}" || exit
-wmake -check-dir "$WM_PROJECT_DIR" 2>/dev/null || {
-    echo "Error (${0##*/}) : not located in \$WM_PROJECT_DIR"
-    echo "    Check your OpenFOAM environment and installation"
-    exit 1
-}
-if [ -f "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments ]
-then  . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments -no-recursion "$@" || \
-    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
-
-#------------------------------------------------------------------------------
-# Additional components
-
-case "$FOAM_MODULE_PREFIX" in
-(false | none)
-    echo ========================================
-    echo "OpenFOAM plugins disabled (prefix=${FOAM_MODULE_PREFIX})"
-    echo "Can be built separately:"
-    echo
-    echo "    ./Allwmake-plugins -prefix=..."
-    echo
-    echo ========================================
-    echo
-    ;;
-(*)
-    # Use wmake -all instead of Allwmake to allow for overrides
-    ( cd "$WM_PROJECT_DIR/plugins" 2>/dev/null && wmake -all )
-esac
-
-#------------------------------------------------------------------------------

CONTRIBUTORS.md

@@ -5,29 +5,22 @@ 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
-- Crist<73>bal Ib<49><62>ez
 - Yoshiaki Inoue
 - Mattijs Janssens
 - Andrew Jackson
@@ -35,7 +28,6 @@ It is likely incomplete...
 - Alexander Kabat vel Job
 - Thilo Knacke
 - Shannon Leakey
-- Sergey Lesnik
 - Tommaso Lucchini
 - Graham Macpherson
 - Alexey Matveichev
@@ -48,7 +40,6 @@ It is likely incomplete...
 - Victor Olesen
 - Evangelos Papoutsis-Kiachagias
 - Juho Peltola
-- Josep Pocurull
 - Johan Roenby
 - Henrik Rusche
 - Bruno Santos
@@ -59,12 +50,8 @@ 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

META-INFO/api-info

@@ -1,2 +1,2 @@
-api=2412
+api=2206
 patch=0

README.md

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

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

applications/solvers/acoustic/acousticFoam/Make/options

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

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;

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

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

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

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

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

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

applications/solvers/basic/potentialFoam/overPotentialFoam/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
 {

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

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

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

applications/solvers/combustion/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_)
     {

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

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

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

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

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

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

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

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

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

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
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

applications/solvers/combustion/chemFoam/YEqn.H

@@ -6,7 +6,7 @@
         solve
         (
             fvm::ddt(rho, Yi) - chemistry.RR(specieI),
-            "Yi"
+            mesh.solver("Yi")
         );
     }
 }

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)

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

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

applications/solvers/combustion/fireFoam/Make/options

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

applications/solvers/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();
 

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

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
         )

applications/solvers/combustion/reactingFoam/EEqn.H

@@ -21,11 +21,6 @@
       + fvOptions(rho, he)
     );
 
-    if (MRF.active())
-    {
-        EEqn += fvc::div(MRF.phi(), p);
-    }
-
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

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

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

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

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

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

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

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

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

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

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

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

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

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

applications/solvers/compressible/rhoPimpleFoam/EEqn.H

@@ -20,11 +20,6 @@
         fvOptions(rho, he)
     );
 
-    if (MRF.active())
-    {
-        EEqn += fvc::div(MRF.phi(), p);
-    }
-
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/compressible/rhoPimpleFoam/Make/options

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

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/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);
 

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

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

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

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

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

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

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

@@ -25,6 +25,17 @@ surfaceScalarField phiHbyA
     fvc::interpolate(rho)*fvc::flux(HbyA)
 );
 
+if (ddtCorr)
+{
+    surfaceScalarField faceMaskOld
+    (
+        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+    );
+
+    phiHbyA +=
+        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf));
+}
+
 fvc::makeRelative(phiHbyA, rho, U);
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
@@ -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;

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

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

applications/solvers/compressible/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())
         {

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

applications/solvers/compressible/rhoSimpleFoam/EEqn.H

@@ -14,11 +14,6 @@
         fvOptions(rho, he)
     );
 
-    if (MRF.active())
-    {
-        EEqn += fvc::div(MRF.phi(), p);
-    }
-
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/compressible/rhoSimpleFoam/Make/options

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

applications/solvers/compressible/sonicFoam/EEqn.H

@@ -9,11 +9,6 @@
         fvOptions(rho, e)
     );
 
-    if (MRF.active())
-    {
-        EEqn += fvc::div(MRF.phi(), p);
-    }
-
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

@@ -127,7 +127,7 @@ int main(int argc, char *argv[])
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
-                    pEqn.solve(p.select(piso.finalInnerIter()));
+                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
                     if (piso.finalNonOrthogonalIter())
                     {
@@ -167,7 +167,7 @@ int main(int argc, char *argv[])
                     fvm::laplacian(rABf, pB) == fvc::div(phiB)
                 );
 
-                pBEqn.solve(pB.select(bpiso.finalInnerIter()));
+                pBEqn.solve(mesh.solver(pB.select(bpiso.finalInnerIter())));
 
                 if (bpiso.finalNonOrthogonalIter())
                 {

applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H

@@ -4,10 +4,10 @@
         sqrt
         (
             2*M_PI*sigma*sqr(aMesh.edgeInterpolation::deltaCoeffs())
-            *mag(aMesh.edgeInterpolation::deltaCoeffs())
+            *aMesh.edgeInterpolation::deltaCoeffs()
             /rhol
         )
-    ).value()*runTime.deltaTValue();
+    ).value()*runTime.deltaT().value();
 
     Info<< "Max Capillary Courant Number = " << CoNumSigma << '\n' << endl;
 }

applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H

@@ -5,7 +5,7 @@
         (
             "h",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
@@ -20,7 +20,7 @@
         (
             "Us",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
@@ -34,7 +34,7 @@
         (
             "phis",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
@@ -48,7 +48,7 @@
         (
             "phi2s",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
@@ -67,7 +67,7 @@
         (
             "Sm",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
@@ -82,7 +82,7 @@
         (
             "Sd",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
@@ -96,7 +96,7 @@
         (
             "Sg",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
@@ -112,7 +112,7 @@
         (
             "ps",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -134,7 +134,7 @@
         (
             "manningField",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
@@ -147,7 +147,7 @@
         (
             "frictionFactor",
             runTime.timeName(),
-            aMesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),

applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H

@@ -1,13 +1,10 @@
-// Volume-to surface mapping object
-const volSurfaceMapping vsm(aMesh);
-
 volVectorField U
 (
     IOobject
     (
         "U",
         runTime.timeName(),
-        mesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
@@ -22,10 +19,13 @@ volScalarField H
     (
         "H",
         runTime.timeName(),
-        mesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
     mesh,
     dimensionedScalar(dimLength, Zero)
 );
+
+// Create volume-to surface mapping object
+volSurfaceMapping vsm(aMesh);

applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C

@@ -87,13 +87,7 @@ int main(int argc, char *argv[])
             (
                 fam::ddt(h, Us)
               + fam::div(phi2s, Us)
-              + fam::Sp
-                (
-                    0.0125
-                   *frictionFactor.internalField()
-                   *mag(Us.internalField()),
-                    Us
-                )
+              + fam::Sp(0.0125*frictionFactor*mag(Us), Us)
              ==
                 Gs*h
               - fam::Sp(Sd, Us)

applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H

@@ -47,10 +47,10 @@ if (aMesh.nInternalEdges())
     );
 
     CoNum = max(SfUfbyDelta/aMesh.magLe())
-        .value()*runTime.deltaTValue();
+        .value()*runTime.deltaT().value();
 
     meanCoNum = (sum(SfUfbyDelta)/sum(aMesh.magLe()))
-        .value()*runTime.deltaTValue();
+        .value()*runTime.deltaT().value();
 
     velMag = max(mag(phis)/aMesh.magLe()).value();
 }

applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H

@@ -5,7 +5,7 @@ areaScalarField Cs
     (
         "Cs",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::MUST_READ,
         IOobject::AUTO_WRITE
     ),
@@ -26,7 +26,7 @@ areaVectorField Us
     (
         "Us",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),
@@ -47,7 +47,6 @@ forAll(Us, faceI)
     Us[faceI].z() =
         Uinf.value()*0.25*R[faceI].x()*R[faceI].z()/sqr(mag(R[faceI]));
 }
-Us.boundaryFieldRef().evaluateCoupled<coupledFaPatch>();
 
 Us -= aMesh.faceAreaNormals()*(aMesh.faceAreaNormals() & Us);
 
@@ -58,7 +57,7 @@ edgeScalarField phis
     (
         "phis",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),

applications/solvers/finiteArea/sphereSurfactantFoam/createFaMesh.H

@@ -0,0 +1,2 @@
+    // Create Finite Area mesh
+    faMesh aMesh(mesh);

applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H

@@ -1,5 +1,5 @@
-    // Volume-to surface mapping object
-    const volSurfaceMapping vsm(aMesh);
+    // Create volume-to surface mapping object
+    volSurfaceMapping vsm(aMesh);
 
     volScalarField Cvf
     (
@@ -7,7 +7,7 @@
         (
             "Cvf",
             runTime.timeName(),
-            mesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -24,7 +24,7 @@
         (
             "U",
             runTime.timeName(),
-            mesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),

applications/solvers/finiteArea/surfactantFoam/createFaFields.H

@@ -5,7 +5,7 @@ areaScalarField Cs
     (
         "Cs",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::MUST_READ,
         IOobject::AUTO_WRITE
     ),
@@ -37,7 +37,7 @@ areaVectorField Us
     (
         "Us",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::MUST_READ,
         IOobject::NO_WRITE
     ),
@@ -51,7 +51,7 @@ edgeScalarField phis
     (
         "phis",
         runTime.timeName(),
-        aMesh.thisDb(),
+        mesh,
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),

applications/solvers/finiteArea/surfactantFoam/createFaMesh.H

@@ -0,0 +1,2 @@
+    // Create Finite Area mesh
+    faMesh aMesh(mesh);

applications/solvers/finiteArea/surfactantFoam/createVolFields.H

@@ -1,5 +1,5 @@
-    // Volume-to surface mapping object
-    const volSurfaceMapping vsm(aMesh);
+    // Create volume-to surface mapping object
+    volSurfaceMapping vsm(aMesh);
 
     volScalarField Cvf
     (
@@ -7,7 +7,7 @@
         (
             "Cvf",
             runTime.timeName(),
-            mesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -24,7 +24,7 @@
         (
             "U",
             runTime.timeName(),
-            mesh.thisDb(),
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -27,7 +27,7 @@
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

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

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/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
         const dictionary& d = mesh.solver
             (

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

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

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

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

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

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

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

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

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -52,7 +52,7 @@ while (pimple.correctNonOrthogonal())
         compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
     );
 
-    p_rghEqn.solve(p_rgh.select(pimple.finalInnerIter()));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {

applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options

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

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

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

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

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

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options

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

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

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

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

@@ -108,7 +108,7 @@ forAll(fluidRegions, i)
                 IOobject::READ_IF_PRESENT,
                 IOobject::NO_WRITE
             ),
-            dimensionedScalar(word::null, dimLength, Zero)
+            dimensionedScalar("hRef", dimLength, Zero) // uses name
         )
     );
 

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

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

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

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

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/Make/files

@@ -1,5 +1,5 @@
 derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C
-solid/solidRegionDiffNo.C
+../solid/solidRegionDiffNo.C
 chtMultiRegionTwoPhaseEulerFoam.C
 
 EXE = $(FOAM_APPBIN)/chtMultiRegionTwoPhaseEulerFoam

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

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

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C

@@ -88,29 +88,26 @@ kappa
 
         case mtLookup:
         {
+            if (mesh.foundObject<volScalarField>(kappaName_))
             {
-                const auto* ptr =
-                    mesh.cfindObject<volScalarField>(kappaName_);
-
-                if (ptr)
-                {
-                    return patch().patchField(*ptr);
-                }
+                return patch().lookupPatchField<volScalarField, scalar>
+                (
+                    kappaName_
+                );
             }
+            else if (mesh.foundObject<volSymmTensorField>(kappaName_))
             {
-                const auto* ptr =
-                    mesh.cfindObject<volSymmTensorField>(kappaName_);
+                const symmTensorField& KWall =
+                    patch().lookupPatchField<volSymmTensorField, scalar>
+                    (
+                        kappaName_
+                    );
 
-                if (ptr)
-                {
-                    const symmTensorField& KWall = patch().patchField(*ptr);
+                const vectorField n(patch().nf());
 
-                    const vectorField n(patch().nf());
-
-                    return n & KWall & n;
-                }
+                return n & KWall & n;
             }
-
+            else
             {
                 FatalErrorInFunction
                     << "Did not find field " << kappaName_
@@ -120,6 +117,9 @@ kappa
                     << " or volSymmTensorField."
                     << exit(FatalError);
             }
+
+
+
             break;
         }
 
@@ -131,8 +131,10 @@ kappa
                 mesh.lookupObject<phaseSystem>("phaseProperties")
             );
 
-            auto tkappaEff = tmp<scalarField>::New(patch().size(), Zero);
-            auto& kappaEff = tkappaEff.ref();
+            tmp<scalarField> kappaEff
+            (
+                new scalarField(patch().size(), 0.0)
+            );
 
             forAll(fluid.phases(), phasei)
             {
@@ -140,10 +142,10 @@ kappa
 
                 const fvPatchScalarField& alpha = phase.boundaryField()[patchi];
 
-                kappaEff += alpha*phase.kappaEff(patchi)();
+                kappaEff.ref() += alpha*phase.kappaEff(patchi)();
             }
 
-            return tkappaEff;
+            return kappaEff;
 
             break;
         }
@@ -159,11 +161,9 @@ kappa
         }
     }
 
-    // Return zero-sized (not nullptr)
-    return tmp<scalarField>::New();
+    return scalarField(0);
 }
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
@@ -243,12 +243,14 @@ turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField
             << exit(FatalError);
     }
 
+    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
 
-    this->readValueEntry(dict, IOobjectOption::MUST_READ);
-
-    if (this->readMixedEntries(dict))
+    if (dict.found("refValue"))
     {
         // Full restart
+        refValue() = scalarField("refValue", dict, p.size());
+        refGrad() = scalarField("refGradient", dict, p.size());
+        valueFraction() = scalarField("valueFraction", dict, p.size());
     }
     else
     {
@@ -292,7 +294,8 @@ updateCoeffs()
 
     // Since we're inside initEvaluate/evaluate there might be processor
     // comms underway. Change the tag we use.
-    const int oldTag = UPstream::incrMsgType();
+    int oldTag = UPstream::msgType();
+    UPstream::msgType() = oldTag+1;
 
     // Get the coupling information from the mappedPatchBase
     const label patchi = patch().index();
@@ -305,11 +308,12 @@ updateCoeffs()
 
     scalarField& Tp = *this;
 
-    const auto& nbrField =
-        refCast
-        <
-            const turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField
-        >(nbrPatch.lookupPatchField<volScalarField>(TnbrName_));
+    const turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField&
+        nbrField = refCast
+            <const turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField>
+            (
+                nbrPatch.lookupPatchField<volScalarField, scalar>(TnbrName_)
+            );
 
     // Swap to obtain full local values of neighbour internal field
     scalarField TcNbr(nbrField.patchInternalField());
@@ -326,13 +330,13 @@ updateCoeffs()
     scalarField qr(Tp.size(), 0.0);
     if (qrName_ != "none")
     {
-        qr = patch().lookupPatchField<volScalarField>(qrName_);
+        qr = patch().lookupPatchField<volScalarField, scalar>(qrName_);
     }
 
     scalarField qrNbr(Tp.size(), 0.0);
     if (qrNbrName_ != "none")
     {
-        qrNbr = nbrPatch.lookupPatchField<volScalarField>(qrNbrName_);
+        qrNbr = nbrPatch.lookupPatchField<volScalarField, scalar>(qrNbrName_);
         mpp.distribute(qrNbr);
     }
 
@@ -470,9 +474,10 @@ updateCoeffs()
             << regionTypeNames_ << nl << exit(FatalError);
     }
 
-    UPstream::msgType(oldTag);  // Restore tag
-
     mixedFvPatchScalarField::updateCoeffs();
+
+    // Restore tag
+    UPstream::msgType() = oldTag;
 }
 
 
@@ -481,7 +486,7 @@ void turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField::write
     Ostream& os
 ) const
 {
-    mixedFvPatchField<scalar>::write(os);
+    mixedFvPatchScalarField::write(os);
     os.writeEntry("kappaMethod", KMethodTypeNames_[method_]);
     os.writeEntryIfDifferent<word>("kappa","none", kappaName_);
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.H

@@ -196,6 +196,18 @@ public:
             const fvPatchFieldMapper&
         );
 
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField
+                (
+                    *this
+                )
+            );
+        }
+
         //- Construct as copy setting internal field reference
         turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField
         (
@@ -203,19 +215,20 @@ 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 turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
         }
 
 

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

@@ -221,7 +221,7 @@ forAll(fluidRegions, i)
                 IOobject::READ_IF_PRESENT,
                 IOobject::NO_WRITE
             ),
-            dimensionedScalar(word::null, dimLength, Zero)
+            dimensionedScalar("hRef", dimLength, Zero)
         )
     );
 

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

@@ -13,11 +13,11 @@ forAll(cumulativeContErrIO, i)
                 "cumulativeContErr",
                 runTime.timeName(),
                 "uniform",
-                mesh.thisDb(),
+                mesh,
                 IOobject::READ_IF_PRESENT,
                 IOobject::AUTO_WRITE
             ),
-            dimensionedScalar(word::null, dimless, Zero)
+            dimensionedScalar(dimless, Zero)
         )
     );
 }

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

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

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

@@ -1,3 +1,5 @@
+    fvMesh& mesh = fluidRegions[i];
+
     twoPhaseSystem& fluid  = phaseSystemFluid[i];
 
     phaseModel& phase1 = fluid.phase1();
@@ -8,13 +10,11 @@
 
     volVectorField& U1 = phase1.URef();
     surfaceScalarField& phi1 = phase1.phiRef();
-    const tmp<surfaceScalarField> talphaPhi1 = phase1.alphaPhi();
-    const auto& alphaPhi1 = talphaPhi1();
+    const surfaceScalarField& alphaPhi1 = phase1.alphaPhi();
 
     volVectorField& U2 = phase2.URef();
     surfaceScalarField& phi2 = phase2.phiRef();
-    const tmp<surfaceScalarField> talphaPhi2 = phase2.alphaPhi();
-    const auto& alphaPhi2 = talphaPhi2();
+    const surfaceScalarField& alphaPhi2 = phase2.alphaPhi();
 
     surfaceScalarField& phi = fluid.phi();