CONFIG: bump patch level (240625)

- the old logic relied on the presence/absence of a comma to decide
  whether to parse as fixed or free format. This logic is faulty when
  handling the final (trailing) entry and will generally lead to the
  last field being truncated when read in.
  Now the caller decides on fixed vs free.

FIX: inconsistent Nastran surface output format

- use FREE format by default. Previously had an odd mix of SHORT
  format when created without options and LONG format (as default)
  when created with format options.

.gitlab/issue_templates/bug.md

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

CONTRIBUTORS.md

@@ -5,8 +5,6 @@ It is likely incomplete...
 
 ## Contributors (alphabetical by surname)
 
-- Tetsuo Aoyagi
-- Akira Azami
 - William Bainbridge
 - Gabriel Barajas
 - Kutalmis Bercin
@@ -21,7 +19,6 @@ It is likely incomplete...
 - Bernhard Gschaider
 - Andrew Heather
 - David Hill
-- Yoshiaki Inoue
 - Mattijs Janssens
 - Andrew Jackson
 - Hrvoje Jasak

META-INFO/api-info

@@ -1,2 +1,2 @@
-api=2301
-patch=230110
+api=2206
+patch=240625

README.md

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

applications/solvers/acoustic/acousticFoam/Make/options

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

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

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

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

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

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

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

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

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

@@ -39,13 +39,13 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             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/reactingFoam/YEqn.H

@@ -38,11 +38,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             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/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C

@@ -165,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 =
@@ -203,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_);
@@ -213,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/U/maxwellSlipUFvPatchVectorField.C

@@ -155,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
     (
@@ -184,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));
@@ -197,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_);
@@ -212,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/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/rhoPimpleFoam/Make/options

@@ -26,5 +26,4 @@ EXE_LIBS = \
     -ltopoChangerFvMesh \
     -lsampling \
     -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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

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

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

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

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

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

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

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

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

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

applications/solvers/compressible/rhoSimpleFoam/Make/options

@@ -21,5 +21,4 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModels \
     -lthermoTools \
     -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H

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

applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H

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

applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H

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

applications/solvers/finiteArea/surfactantFoam/createVolFields.H

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

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

@@ -28,5 +28,4 @@ EXE_LIBS = \
     -ldynamicMesh \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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

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

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

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

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

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

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

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

applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options

@@ -22,5 +22,4 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModels \
     -lthermoTools \
     -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -40,5 +40,4 @@ EXE_LIBS = \
     -lthermoTools \
     -lradiationModels \
     -lregionModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

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

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

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

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

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

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

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

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

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

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

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

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/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);
 }
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
@@ -308,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());
@@ -329,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);
     }
 
@@ -485,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/fluid/readFluidMultiRegionPIMPLEControls.H

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

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

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

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

@@ -48,7 +48,7 @@ if (Y.size())
 
             fvOptions.correct(Yi);
 
-            Yi.clamp_min(0);
+            Yi.max(0.0);
             Yt += Yi;
         }
     }
@@ -56,6 +56,6 @@ if (Y.size())
     if (Y.size())
     {
         Y[inertIndex] = scalar(1) - Yt;
-        Y[inertIndex].clamp_min(0);
+        Y[inertIndex].max(0.0);
     }
 }

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

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

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

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

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

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

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

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

applications/solvers/heatTransfer/solidFoam/createFields.H

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

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointOutletPressure/adjointOutletPressureFvPatchScalarField.C

@@ -89,10 +89,17 @@ void Foam::adjointOutletPressureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const auto& phip = patch().lookupPatchField<surfaceScalarField>("phi");
-    const auto& phiap = patch().lookupPatchField<surfaceScalarField>("phia");
-    const auto& Up = patch().lookupPatchField<volVectorField>("U");
-    const auto& Uap = patch().lookupPatchField<volVectorField>("Ua");
+    const fvsPatchField<scalar>& phip =
+        patch().lookupPatchField<surfaceScalarField, scalar>("phi");
+
+    const fvsPatchField<scalar>& phiap =
+        patch().lookupPatchField<surfaceScalarField, scalar>("phia");
+
+    const fvPatchField<vector>& Up =
+        patch().lookupPatchField<volVectorField, vector>("U");
+
+    const fvPatchField<vector>& Uap =
+        patch().lookupPatchField<volVectorField, vector>("Ua");
 
     operator==((phiap/patch().magSf() - 1.0)*phip/patch().magSf() + (Up & Uap));
 
@@ -102,8 +109,8 @@ void Foam::adjointOutletPressureFvPatchScalarField::updateCoeffs()
 
 void Foam::adjointOutletPressureFvPatchScalarField::write(Ostream& os) const
 {
-    fvPatchField<scalar>::write(os);
-    fvPatchField<scalar>::writeValueEntry(os);
+    fvPatchScalarField::write(os);
+    writeEntry("value", os);
 }
 
 

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C

@@ -90,8 +90,11 @@ void Foam::adjointOutletVelocityFvPatchVectorField::updateCoeffs()
         return;
     }
 
-    const auto& phiap = patch().lookupPatchField<surfaceScalarField>("phia");
-    const auto& Up = patch().lookupPatchField<volVectorField>("U");
+    const fvsPatchField<scalar>& phiap =
+        patch().lookupPatchField<surfaceScalarField, scalar>("phia");
+
+    const fvPatchField<vector>& Up =
+        patch().lookupPatchField<volVectorField, vector>("U");
 
     scalarField Un(mag(patch().nf() & Up));
     vectorField UtHat((Up - patch().nf()*Un)/(Un + SMALL));
@@ -107,8 +110,8 @@ void Foam::adjointOutletVelocityFvPatchVectorField::updateCoeffs()
 
 void Foam::adjointOutletVelocityFvPatchVectorField::write(Ostream& os) const
 {
-    fvPatchField<vector>::write(os);
-    fvPatchField<vector>::writeValueEntry(os);
+    fvPatchVectorField::write(os);
+    writeEntry("value", os);
 }
 
 

applications/solvers/incompressible/pimpleFoam/Make/options

@@ -22,5 +22,4 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
     -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels

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

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

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

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

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

@@ -94,7 +94,7 @@ if (mesh.changing())
             {
                 if (refCells[zoneId] != -1)
                 {
-                    validCells.push_back(refCells[zoneId]);
+                    validCells.append(refCells[zoneId]);
                 }
             }
 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

applications/solvers/lagrangian/coalChemistryFoam/Make/options

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

applications/solvers/lagrangian/coalChemistryFoam/YEqn.H

@@ -40,11 +40,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             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/lagrangian/reactingParcelFoam/Make/options

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

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -41,11 +41,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             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/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options

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

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

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

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

@@ -38,11 +38,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             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/lagrangian/simpleCoalParcelFoam/Make/options

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

applications/solvers/lagrangian/simpleCoalParcelFoam/YEqn.H

@@ -38,11 +38,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             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/lagrangian/sprayFoam/Make/options

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

applications/solvers/lagrangian/sprayFoam/YEqn.H

@@ -39,11 +39,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
             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/lagrangian/sprayFoam/engineFoam/Make/options

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

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

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

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

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

applications/solvers/multiphase/MPPICInterFoam/Make/options

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

applications/solvers/multiphase/cavitatingFoam/alphavPsi.H

@@ -1,5 +1,14 @@
 {
-    alphav = clamp((rho - rholSat)/(rhovSat - rholSat), zero_one{});
+    alphav =
+        max
+        (
+            min
+            (
+                (rho - rholSat)/(rhovSat - rholSat),
+                scalar(1)
+            ),
+            scalar(0)
+        );
     alphal = 1.0 - alphav;
 
     Info<< "max-min alphav: " << max(alphav).value()

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

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

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

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

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

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

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C

@@ -126,9 +126,9 @@ alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
 
 void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
 {
-    fvPatchField<scalar>::write(os);
+    fvPatchScalarField::write(os);
     os.writeEntry("thetaProperties", thetaProps_);
-    fvPatchField<scalar>::writeValueEntry(os);
+    writeEntry("value", os);
 }
 
 

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C

@@ -779,7 +779,7 @@ Foam::multiphaseMixtureThermo::surfaceTensionForce() const
 
             auto sigma = sigmas_.cfind(interfacePair(alpha1, alpha2));
 
-            if (!sigma.good())
+            if (!sigma.found())
             {
                 FatalErrorInFunction
                     << "Cannot find interface " << interfacePair(alpha1, alpha2)
@@ -907,7 +907,7 @@ void Foam::multiphaseMixtureThermo::correctContactAngle
             const auto tp =
                 acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
 
-            if (!tp.good())
+            if (!tp.found())
             {
                 FatalErrorInFunction
                     << "Cannot find interface " << interfacePair(alpha1, alpha2)

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.H

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_DTRMParticle_H
-#define Foam_DTRMParticle_H
+#ifndef DTRMParticle_H
+#define DTRMParticle_H
 
 #include "particle.H"
 #include "IOstream.H"
@@ -50,13 +50,12 @@ SourceFiles
 namespace Foam
 {
 
-// Forward Declarations
 class DTRMParticle;
-
-Ostream& operator<<(Ostream&, const DTRMParticle&);
-
 using namespace Foam::radiation;
 
+// Forward declaration of friend functions
+Ostream& operator<<(Ostream&, const DTRMParticle&);
+
 /*---------------------------------------------------------------------------*\
                         Class DTRMParticle Declaration
 \*---------------------------------------------------------------------------*/
@@ -131,7 +130,7 @@ public:
             );
 
 
-        // Member Functions
+        // Member functions
 
             inline const interpolationCell<scalar>& aInterp() const;
             inline const interpolationCell<scalar>& eInterp() const;
@@ -233,34 +232,37 @@ public:
     // Access
 
         //- Return const access to the initial position
-        const point& p0() const noexcept { return p0_; }
+        inline const point& p0() const;
 
         //- Return const access to the target position
-        const point& p1() const noexcept { return p1_; }
+        inline const point& p1() const;
 
         //- Return const access to the initial intensity
-        scalar I0() const noexcept { return I0_; }
+        inline scalar I0() const;
 
         //- Return const access to the current intensity
-        scalar I() const noexcept { return I_; }
+        inline scalar I() const;
 
         //- Return const access dA
-        scalar dA() const noexcept { return dA_; }
+        inline scalar dA() const;
 
 
     // Edit
 
         //- Return access to the target position
-        point& p1() noexcept { return p1_; }
+        inline point& p1();
 
         //- Return access to the initial intensity
-        scalar& I0() noexcept { return I0_; }
+        inline scalar& I0();
 
         //- Return access to the current intensity
-        scalar& I() noexcept { return I_; }
+        inline scalar& I();
 
         //- Return access to dA
-        scalar& dA() noexcept { return dA_; }
+        inline scalar& dA();
+
+        //- Return access to reflectedId
+        inline label& reflectedId();
 
 
    // Tracking

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleI.H

@@ -107,4 +107,58 @@ inline Foam::scalar& Foam::DTRMParticle::trackingData::Q(label celli)
 }
 
 
+inline const Foam::point& Foam::DTRMParticle::p0() const
+{
+    return p0_;
+}
+
+
+inline const Foam::point& Foam::DTRMParticle::p1() const
+{
+    return p1_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::I0() const
+{
+    return I0_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::I() const
+{
+    return I_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::dA() const
+{
+    return dA_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::dA()
+{
+    return dA_;
+}
+
+
+inline Foam::point& Foam::DTRMParticle::p1()
+{
+    return p1_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::I0()
+{
+    return I0_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::I()
+{
+    return I_;
+}
+
+
 // ************************************************************************* //

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

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

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

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

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.C

@@ -99,6 +99,17 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotAlphal() const
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
 {
+
+    volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
+
+    volScalarField limitedAlpha2
+    (
+        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
+    );
+
     const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
 
     const twoPhaseMixtureEThermo& thermo =
@@ -113,15 +124,11 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
 
     volScalarField mDotE
     (
-        "mDotE",
-        coeffE_*mixture_.rho1()*clamp(mixture_.alpha1(), zero_one{})
-      * max(T - TSat, T0)
+        "mDotE", coeffE_*mixture_.rho1()*limitedAlpha1*max(T - TSat, T0)
     );
     volScalarField mDotC
     (
-        "mDotC",
-        coeffC_*mixture_.rho2()*clamp(mixture_.alpha2(), zero_one{})
-      * max(TSat - T, T0)
+        "mDotC", coeffC_*mixture_.rho2()*limitedAlpha2*max(TSat - T, T0)
     );
 
     if (mesh_.time().outputTime())
@@ -141,6 +148,16 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotDeltaT() const
 {
+    volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
+
+    volScalarField limitedAlpha2
+    (
+        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
+    );
+
     const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
 
     const twoPhaseMixtureEThermo& thermo =
@@ -153,14 +170,8 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotDeltaT() const
 
     return Pair<tmp<volScalarField>>
     (
-        (
-            coeffC_*mixture_.rho2()*clamp(mixture_.alpha2(), zero_one{})
-          * pos(TSat - T)
-        ),
-        (
-            coeffE_*mixture_.rho1()*clamp(mixture_.alpha1(), zero_one{})
-          * pos(T - TSat)
-        )
+        coeffC_*mixture_.rho2()*limitedAlpha2*pos(TSat - T),
+        coeffE_*mixture_.rho1()*limitedAlpha1*pos(T - TSat)
     );
 }
 
@@ -190,17 +201,25 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
 
     const dimensionedScalar& TSat = thermo.TSat();
 
-    const dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
+    dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
+
+    volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
+
+    volScalarField limitedAlpha2
+    (
+        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
+    );
 
     const volScalarField Vcoeff
     (
-        coeffE_*mixture_.rho1()*clamp(mixture_.alpha1(), zero_one{})
-      * L*pos(T - TSat)
+        coeffE_*mixture_.rho1()*limitedAlpha1*L*pos(T - TSat)
     );
     const volScalarField Ccoeff
     (
-        coeffC_*mixture_.rho2()*clamp(mixture_.alpha2(), zero_one{})
-      * L*pos(TSat - T)
+        coeffC_*mixture_.rho2()*limitedAlpha2*L*pos(TSat - T)
     );
 
     TSource =

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C

@@ -167,10 +167,20 @@ Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
 mDotAlphal() const
 {
+    volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
+
+    volScalarField limitedAlpha2
+    (
+        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
+    );
+
     return Pair<tmp<volScalarField>>
     (
-        (mDotc_/clamp(mixture_.alpha2(), scalarMinMax(SMALL, 1))),
-       -(mDote_/clamp(mixture_.alpha1(), scalarMinMax(SMALL, 1)))
+        (mDotc_/(limitedAlpha2 + SMALL)),
+       -(mDote_/(limitedAlpha1 + SMALL))
     );
 }
 

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C

@@ -155,7 +155,7 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::Cp() const
 {
     const volScalarField limitedAlpha1
     (
-        clamp(alpha1_, zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
     return tmp<volScalarField>
@@ -176,11 +176,13 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::Cp
     const label patchi
 ) const
 {
-    const scalarField alpha1p
+    const volScalarField limitedAlpha1
     (
-        clamp(alpha1_.boundaryField()[patchi], zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
+    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
+
     return
     (
         alpha1p*Cp1().value() + (scalar(1) - alpha1p)*Cp2().value()
@@ -192,7 +194,7 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::rho() const
 {
     const volScalarField limitedAlpha1
     (
-        clamp(alpha1_, zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
     return tmp<volScalarField>
@@ -212,11 +214,13 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::rho
     const label patchi
 ) const
 {
-    const scalarField alpha1p
+    const volScalarField limitedAlpha1
     (
-        clamp(alpha1_.boundaryField()[patchi], zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
+    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
+
     return
     (
         alpha1p*rho1().value() + (scalar(1) - alpha1p)*rho2().value()
@@ -228,7 +232,7 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::Cv() const
 {
     const volScalarField limitedAlpha1
     (
-        clamp(alpha1_, zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
     return tmp<volScalarField>
@@ -249,11 +253,13 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::Cv
     const label patchi
 ) const
 {
-    const scalarField alpha1p
+    const volScalarField limitedAlpha1
     (
-        clamp(alpha1_.boundaryField()[patchi], zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
+    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
+
     return
     (
         alpha1p*Cv1().value() + (scalar(1) - alpha1p)*Cv2().value()
@@ -333,7 +339,7 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::kappa() const
 {
     const volScalarField limitedAlpha1
     (
-        clamp(alpha1_, zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
     return tmp<volScalarField>
@@ -352,11 +358,13 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::kappa
     const label patchi
 ) const
 {
-    const scalarField alpha1p
+    const volScalarField limitedAlpha1
     (
-        clamp(alpha1_.boundaryField()[patchi], zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
+    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
+
     return (alpha1p*kappa1().value() + (1 - alpha1p)*kappa2().value());
 }
 
@@ -394,11 +402,13 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::kappaEff
     const label patchi
 ) const
 {
-    const scalarField alpha1p
+    const volScalarField limitedAlpha1
     (
-        clamp(alpha1_.boundaryField()[patchi], zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
+    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
+
     return
         (alpha1p*kappa1().value() + (1 - alpha1p)*kappa2().value()) + kappat;
 
@@ -425,11 +435,13 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::alphaEff
     const label patchi
 ) const
 {
-    const scalarField alpha1p
+    const volScalarField limitedAlpha1
     (
-        clamp(alpha1_.boundaryField()[patchi], zero_one{})
+        min(max(alpha1_, scalar(0)), scalar(1))
     );
 
+    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
+
     const scalarField rho
     (
         alpha1p*rho1().value() + (1.0 - alpha1p)*rho2().value()

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2011 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,10 +36,7 @@ scalar meanCoNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField phiMask
-    (
-        localMin<scalar>(mesh).interpolate(cellMask + interpolatedCells)
-    );
+    surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     scalarField sumPhi
     (

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

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

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

@@ -103,7 +103,7 @@
             {
                 if (refCells[zoneId] != -1)
                 {
-                    validCells.push_back(refCells[zoneId]);
+                    validCells.append(refCells[zoneId]);
                 }
             }
 

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

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

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

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

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

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

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

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

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

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