Feature evaluation check

etc/controlDict

@@ -213,6 +213,10 @@ OptimisationSwitches
     //  in commit da787200.  Default is to use the formulation from v1712
     //  see ddtScheme.C
     experimentalDdtCorr 0;
+
+    //- Enable enforced consistency of constraint bcs after 'local' operations.
+    //  Default is on. Set to 0/false to revert to <v2306 behaviour
+    //localConsistency 0;
 }
 
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017,2022 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,6 +33,159 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<class CheckPatchFieldType>
+bool Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::checkConsistency
+(
+    const scalar tol,
+    const bool doExit
+) const
+{
+    if (!this->size())
+    {
+        return true;
+    }
+
+    if (debug&2)
+    {
+        const auto& pfld0 = this->operator[](0);
+        PoutInFunction
+            << " Checking boundary consistency for field "
+            << pfld0.internalField().name()
+            << endl;
+    }
+
+    auto& bfld = const_cast<GeometricBoundaryField<Type, PatchField, GeoMesh>&>
+    (
+        *this
+    );
+
+
+    // Store old value
+    List<Field<Type>> oldBfld(this->size());
+    boolList oldUpdated(this->size());
+    //Note: areaFields (finiteArea) do not have manipulatedMatrix() flag. TBD.
+    //boolList oldManipulated(this->size());
+
+    for (auto& pfld : bfld)
+    {
+        if (isA<CheckPatchFieldType>(pfld))
+        {
+            const label patchi = pfld.patch().index();
+            oldUpdated[patchi] = pfld.updated();
+            oldBfld[patchi] = pfld;
+            //oldManipulated[patchi] = pfld.manipulatedMatrix();
+        }
+    }
+
+
+    // Re-evaluate
+    {
+        const label startOfRequests = UPstream::nRequests();
+
+        for (auto& pfld : bfld)
+        {
+            if (isA<CheckPatchFieldType>(pfld))
+            {
+                pfld.initEvaluate(UPstream::commsTypes::nonBlocking);
+            }
+        }
+
+        // Wait for outstanding requests
+        UPstream::waitRequests(startOfRequests);
+
+        for (auto& pfld : bfld)
+        {
+            if (isA<CheckPatchFieldType>(pfld))
+            {
+                pfld.evaluate(UPstream::commsTypes::nonBlocking);
+            }
+        }
+    }
+
+
+    // Check
+    bool ok = true;
+    for (auto& pfld : bfld)
+    {
+        if (isA<CheckPatchFieldType>(pfld))
+        {
+            const label patchi = pfld.patch().index();
+            const auto& oldPfld = oldBfld[patchi];
+
+            forAll(pfld, facei)
+            {
+                if (mag(pfld[facei]-oldPfld[facei]) > tol)
+                {
+                    ok = false;
+                    break;
+                }
+            }
+
+            if (!ok)
+            {
+                if (doExit)
+                {
+                    FatalErrorInFunction << "Field "
+                        << pfld.internalField().name()
+                        << " is not evaluated?"
+                        << " On patch " << pfld.patch().name()
+                        << " type " << pfld.type()
+                        << " : average of field = "
+                        << average(oldPfld)
+                        << ". Average of evaluated field = "
+                        << average(pfld)
+                        << ". Difference:" << average(pfld-oldPfld)
+                        << ". Tolerance:" << tol
+                        << exit(FatalError);
+                }
+                else
+                {
+                    WarningInFunction << "Field "
+                        << pfld.internalField().name()
+                        << " is not evaluated?"
+                        << " On patch " << pfld.patch().name()
+                        << " type " << pfld.type()
+                        << " : average of field = "
+                        << average(oldPfld)
+                        << ". Average of evaluated field = "
+                        << average(pfld)
+                        << ". Difference:" << average(pfld-oldPfld)
+                        << ". Tolerance:" << tol
+                        << endl;
+
+                    // Skip other patches
+                    break;
+                }
+            }
+        }
+    }
+
+    // Restore bfld, updated
+    for (auto& pfld : bfld)
+    {
+        if (isA<CheckPatchFieldType>(pfld))
+        {
+            const label patchi = pfld.patch().index();
+            pfld.setUpdated(oldUpdated[patchi]);
+            Field<Type>& vals = pfld;
+            vals = std::move(oldBfld[patchi]);
+            //pfld.setManipulated(oldManipulated[patchi]);
+        }
+    }
+
+    if (debug&2)
+    {
+        const auto& pfld0 = this->operator[](0);
+        PoutInFunction
+            << " Result of checking for field "
+            << pfld0.internalField().name() << " : " << ok << endl;
+    }
+
+    return ok;
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::readField
 (
@@ -492,9 +645,79 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluate()
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluateLocal()
+{
+    ///if (GeometricField<Type, PatchField, GeoMesh::debug)
+    ///{
+    ///    InfoInFunction << nl;
+    ///}
+
+    if (!localConsistency)
+    {
+        return;
+    }
+
+    const UPstream::commsTypes commsType = UPstream::defaultCommsType;
+    const label startOfRequests = UPstream::nRequests();
+
+    if
+    (
+        commsType == UPstream::commsTypes::blocking
+     || commsType == UPstream::commsTypes::nonBlocking
+    )
+    {
+        for (auto& pfld : *this)
+        {
+            pfld.initEvaluateLocal(commsType);
+        }
+
+        // Wait for outstanding requests
+        if (commsType == Pstream::commsTypes::nonBlocking)
+        {
+            UPstream::waitRequests(startOfRequests);
+        }
+
+        for (auto& pfld : *this)
+        {
+            pfld.evaluateLocal(commsType);
+        }
+    }
+    else if (commsType == UPstream::commsTypes::scheduled)
+    {
+        const lduSchedule& patchSchedule =
+            bmesh_.mesh().globalData().patchSchedule();
+
+        for (const auto& schedEval : patchSchedule)
+        {
+            const label patchi = schedEval.patch;
+            auto& pfld = (*this)[patchi];
+
+            if (schedEval.init)
+            {
+                pfld.initEvaluateLocal(commsType);
+            }
+            else
+            {
+                pfld.evaluateLocal(commsType);
+            }
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unsupported communications type "
+            << UPstream::commsTypeNames[commsType]
+            << exit(FatalError);
+    }
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 template<class CoupledPatchType>
-void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::evaluateCoupled()
+void
+Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>
+::evaluateCoupled()
 {
     ///if (GeometricField<Type, PatchField, GeoMesh::debug)
     ///{
@@ -677,6 +900,17 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::writeEntries
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+bool Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::check
+(
+) const
+{
+    // Dummy op - template specialisations provide logic (usually call
+    // to checkConsistency)
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Type, template<class> class PatchField, class GeoMesh>

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017,2022 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,8 +86,27 @@ private:
         const BoundaryMesh& bmesh_;
 
 
+    // Private Member Functions
+
+        //- Helper: check if patchfields have been evaluated. If not:
+        //      exit = true  : FatalError
+        //      exit = false : return bool
+        template<class CheckPatchField>
+        bool checkConsistency(const scalar tol, const bool exitIfBad) const;
+
+
 public:
 
+    //- Enable debug
+    static int debug;
+
+    //- User-defined tolerance (for consistency checks)
+    static scalar tolerance;
+
+    //- Enable local consistency
+    static int localConsistency;
+
+
     // Constructors
 
         //- Construct from a BoundaryMesh, setting patches later
@@ -170,6 +189,9 @@ public:
         //- Evaluate boundary conditions
         void evaluate();
 
+        //- Evaluate boundary conditions after change in local values
+        void evaluateLocal();
+
         //- Evaluate boundary conditions on a subset of coupled patches
         template<class CoupledPatchType>
         void evaluateCoupled();
@@ -194,6 +216,9 @@ public:
         //- Write dictionary entries of the individual boundary fields.
         void writeEntries(Ostream& os) const;
 
+        //- Helper: check if field has been evaluated. See instantiations.
+        bool check() const;
+
 
     // Member Operators
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -887,6 +887,16 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::clone() const
 template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::GeometricField<Type, PatchField, GeoMesh>::~GeometricField()
 {
+    /*
+    if (debug)
+    {
+        // Problem: temporary fields might have their internal field
+        // already stolen so boundary fields will not be able to access the
+        // internal field anymore
+        boundaryField_.check();
+    }
+    */
+
     deleteDemandDrivenData(field0Ptr_);
     deleteDemandDrivenData(fieldPrevIterPtr_);
 }
@@ -1097,6 +1107,16 @@ correctBoundaryConditions()
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::
+correctLocalBoundaryConditions()
+{
+    this->setUpToDate();
+    storeOldTimes();
+    boundaryField_.evaluateLocal();
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 bool Foam::GeometricField<Type, PatchField, GeoMesh>::needReference() const
 {

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H

@@ -635,6 +635,10 @@ public:
         //- Correct boundary field
         void correctBoundaryConditions();
 
+        //- Correct boundary conditions after a purely local operation. Is
+        //  dummy for e.g. processor boundary conditions
+        void correctLocalBoundaryConditions();
+
         //- Does the field need a reference level for solution
         bool needReference() const;
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C

@@ -55,6 +55,10 @@ void component
     component(result.primitiveFieldRef(), f1.primitiveField(), d);
     component(result.boundaryFieldRef(), f1.boundaryField(), d);
     result.oriented() = f1.oriented();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -68,6 +72,10 @@ void T
     T(result.primitiveFieldRef(), f1.primitiveField());
     T(result.boundaryFieldRef(), f1.boundaryField());
     result.oriented() = f1.oriented();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -88,6 +96,11 @@ void pow
     pow(result.primitiveFieldRef(), f1.primitiveField(), r);
     pow(result.boundaryFieldRef(), f1.boundaryField(), r);
     result.oriented() = pow(f1.oriented(), r);
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -167,6 +180,11 @@ void sqr
     sqr(result.primitiveFieldRef(), f1.primitiveField());
     sqr(result.boundaryFieldRef(), f1.boundaryField());
     result.oriented() = sqr(f1.oriented());
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -239,6 +257,11 @@ void magSqr
     magSqr(result.primitiveFieldRef(), f1.primitiveField());
     magSqr(result.boundaryFieldRef(), f1.boundaryField());
     result.oriented() = magSqr(f1.oriented());
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -288,6 +311,11 @@ void mag
     mag(result.primitiveFieldRef(), f1.primitiveField());
     mag(result.boundaryFieldRef(), f1.boundaryField());
     result.oriented() = mag(f1.oriented());
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -342,6 +370,10 @@ void cmptAv
     cmptAv(result.primitiveFieldRef(), f1.primitiveField());
     cmptAv(result.boundaryFieldRef(), f1.boundaryField());
     result.oriented() = cmptAv(f1.oriented());
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -497,6 +529,11 @@ void clamp
     clamp(result.primitiveFieldRef(), f1.primitiveField(), range);
     clamp(result.boundaryFieldRef(), f1.boundaryField(), range);
     result.oriented() = f1.oriented();
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -596,6 +633,11 @@ void OpFunc                                                                    \
     );                                                                         \
                                                                                \
     result.oriented() = (f1.oriented() Op f2.oriented());                      \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<Type1, PatchField, GeoMesh>::debug)             \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
                                                                                \
@@ -734,6 +776,10 @@ void OpFunc                                                                    \
     Foam::OpFunc(result.primitiveFieldRef(), f1.primitiveField(), dvs.value());\
     Foam::OpFunc(result.boundaryFieldRef(), f1.boundaryField(), dvs.value());  \
     result.oriented() = f1.oriented();                                         \
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)              \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
 template                                                                       \
@@ -838,6 +884,10 @@ void OpFunc                                                                    \
     Foam::OpFunc(result.primitiveFieldRef(), dvs.value(), f2.primitiveField());\
     Foam::OpFunc(result.boundaryFieldRef(), dvs.value(), f2.boundaryField());  \
     result.oriented() = f2.oriented();                                         \
+    if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)              \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
 template                                                                       \

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctionsM.C

@@ -47,6 +47,11 @@ void Func                                                                      \
     Foam::Func(result.primitiveFieldRef(), f1.primitiveField());               \
     Foam::Func(result.boundaryFieldRef(), f1.boundaryField());                 \
     result.oriented() = f1.oriented();                                         \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
                                                                                \
@@ -105,6 +110,11 @@ void OpFunc                                                                    \
     Foam::OpFunc(result.primitiveFieldRef(), f1.primitiveField());             \
     Foam::OpFunc(result.boundaryFieldRef(), f1.boundaryField());               \
     result.oriented() = f1.oriented();                                         \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
                                                                                \
@@ -174,6 +184,11 @@ void Func                                                                      \
         f2.boundaryField()                                                     \
     );                                                                         \
     result.oriented() = Func(f1.oriented(), f2.oriented());                    \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
                                                                                \
@@ -286,6 +301,11 @@ void Func                                                                      \
     Foam::Func(result.primitiveFieldRef(), dt1.value(), f2.primitiveField());  \
     Foam::Func(result.boundaryFieldRef(), dt1.value(), f2.boundaryField());    \
     result.oriented() = f2.oriented();                                         \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
                                                                                \
@@ -367,6 +387,11 @@ void Func                                                                      \
     Foam::Func(result.primitiveFieldRef(), f1.primitiveField(), dt2.value());  \
     Foam::Func(result.boundaryFieldRef(), f1.boundaryField(), dt2.value());    \
     result.oriented() = f1.oriented();                                         \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
                                                                                \
@@ -465,6 +490,11 @@ void OpFunc                                                                    \
         f2.boundaryField()                                                     \
     );                                                                         \
     result.oriented() = (f1.oriented() Op f2.oriented());                      \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
                                                                                \
@@ -576,7 +606,11 @@ void OpFunc                                                                    \
     Foam::OpFunc(result.primitiveFieldRef(), dt1.value(), f2.primitiveField());\
     Foam::OpFunc(result.boundaryFieldRef(), dt1.value(), f2.boundaryField());  \
     result.oriented() = f2.oriented();                                         \
-                                                                               \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -657,6 +691,11 @@ void OpFunc                                                                    \
     Foam::OpFunc(result.primitiveFieldRef(), f1.primitiveField(), dt2.value());\
     Foam::OpFunc(result.boundaryFieldRef(), f1.boundaryField(), dt2.value());  \
     result.oriented() = f1.oriented();                                         \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
                                                                                \
@@ -758,6 +797,11 @@ void Func                                                                      \
         f3.boundaryField()                                                     \
     );                                                                         \
     result.oriented() = Func(f1.oriented(), f2.oriented());                    \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \
@@ -990,6 +1034,11 @@ void Func                                                                      \
         dt3.value()                                                            \
     );                                                                         \
     result.oriented() = Func(f1.oriented(), f2.oriented());                    \
+    result.correctLocalBoundaryConditions();                                   \
+    if (GeometricBoundaryField<ReturnType, PatchField, GeoMesh>::debug)        \
+    {                                                                          \
+        result.boundaryField().check();                                        \
+    }                                                                          \
 }                                                                              \
                                                                                \
 TEMPLATE                                                                       \

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldOps.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,6 +86,11 @@ void assign
             op
         );
     }
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Tout, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -127,6 +132,11 @@ void assign
             bop
         );
     }
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<Tout, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -167,6 +177,11 @@ void ternary
             bop
         );
     }
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<T, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 
@@ -213,6 +228,11 @@ void ternarySelect
             flip
         );
     }
+    result.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<T, PatchField, GeoMesh>::debug)
+    {
+        result.boundaryField().check();
+    }
 }
 
 

src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C

@@ -121,6 +121,11 @@ void pow
 {
     pow(Pow.primitiveFieldRef(), gsf1.primitiveField(), gsf2.primitiveField());
     pow(Pow.boundaryFieldRef(), gsf1.boundaryField(), gsf2.boundaryField());
+    Pow.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<scalar, PatchField, GeoMesh>::debug)
+    {
+        Pow.boundaryField().check();
+    }
 }
 
 
@@ -326,6 +331,11 @@ void pow
 {
     pow(tPow.primitiveFieldRef(), gsf.primitiveField(), ds.value());
     pow(tPow.boundaryFieldRef(), gsf.boundaryField(), ds.value());
+    tPow.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<scalar, PatchField, GeoMesh>::debug)
+    {
+        tPow.boundaryField().check();
+    }
 }
 
 
@@ -429,6 +439,11 @@ void pow
 {
     pow(tPow.primitiveFieldRef(), ds.value(), gsf.primitiveField());
     pow(tPow.boundaryFieldRef(), ds.value(), gsf.boundaryField());
+    tPow.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<scalar, PatchField, GeoMesh>::debug)
+    {
+        tPow.boundaryField().check();
+    }
 }
 
 
@@ -565,6 +580,13 @@ void atan2
         gsf1.boundaryField(),
         gsf2.boundaryField()
     );
+
+    Atan2.correctLocalBoundaryConditions();
+
+    if (GeometricBoundaryField<scalar, PatchField, GeoMesh>::debug)
+    {
+        Atan2.boundaryField().check();
+    }
 }
 
 
@@ -682,6 +704,11 @@ void atan2
 {
     atan2(tAtan2.primitiveFieldRef(), gsf.primitiveField(), ds.value());
     atan2(tAtan2.boundaryFieldRef(), gsf.boundaryField(), ds.value());
+    tAtan2.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<scalar, PatchField, GeoMesh>::debug)
+    {
+        tAtan2.boundaryField().check();
+    }
 }
 
 
@@ -761,6 +788,11 @@ void atan2
 {
     atan2(tAtan2.primitiveFieldRef(), ds.value(), gsf.primitiveField());
     atan2(tAtan2.boundaryFieldRef(), ds.value(), gsf.boundaryField());
+    tAtan2.correctLocalBoundaryConditions();
+    if (GeometricBoundaryField<scalar, PatchField, GeoMesh>::debug)
+    {
+        tAtan2.boundaryField().check();
+    }
 }
 
 
@@ -886,6 +918,11 @@ void func                                                                      \
 {                                                                              \
     func(gsf.primitiveFieldRef(), n, gsf1.primitiveField());                   \
     func(gsf.boundaryFieldRef(), n, gsf1.boundaryField());                     \
+    gsf.correctLocalBoundaryConditions();                                      \
+    if (GeometricBoundaryField<scalar, PatchField, GeoMesh>::debug)            \
+    {                                                                          \
+        gsf.boundaryField().check();                                           \
+    }                                                                          \
 }                                                                              \
                                                                                \
 template<template<class> class PatchField, class GeoMesh>                      \

src/OpenFOAM/fields/GeometricFields/pointFields/pointFields.C

@@ -27,6 +27,7 @@ License
 
 #include "polyMesh.H"
 #include "pointFields.H"
+#include "registerSwitch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -45,6 +46,98 @@ defineTemplateTypeNameAndDebug(pointSphericalTensorField, 0);
 defineTemplateTypeNameAndDebug(pointSymmTensorField, 0);
 defineTemplateTypeNameAndDebug(pointTensorField, 0);
 
+
+defineTemplateDebugSwitchWithName
+(
+    pointScalarField::Boundary,
+    "pointScalarField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    pointVectorField::Boundary,
+    "pointVectorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    pointSphericalTensorField::Boundary,
+    "pointSphericalTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    pointSymmTensorField::Boundary,
+    "pointSymmTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    pointTensorField::Boundary,
+    "pointTensorField::Boundary",
+    0
+);
+
+
+
+// Local-ops consistency enforcing
+
+template<> int pointScalarField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "localConsistency",
+    int,
+    Foam::pointScalarField::Boundary::localConsistency
+);
+
+template<> int pointVectorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "pointVectorField::Boundary::localConsistency",
+    int,
+    Foam::pointVectorField::Boundary::localConsistency
+);
+
+template<> int pointSphericalTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "pointSphericalTensorField::Boundary::localConsistency",
+    int,
+    Foam::pointSphericalTensorField::Boundary::localConsistency
+);
+
+template<> int pointSymmTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "pointSymmTensorField::Boundary::localConsistency",
+    int,
+    Foam::pointSymmTensorField::Boundary::localConsistency
+);
+
+template<> int pointTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "pointTensorField::Boundary::localConsistency",
+    int,
+    Foam::pointTensorField::Boundary::localConsistency
+);
+
+
 } // End namespace Foam
 
 

src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H

@@ -527,6 +527,23 @@ public:
             const Pstream::commsTypes commsType = Pstream::commsTypes::blocking
         );
 
+            //- Initialise the evaluation of the patch field after a local
+            //  operation
+            virtual void initEvaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
+
+            //- Evaluate the patch field after a local operation (e.g. *=)
+            virtual void evaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
+
 
     // Other
 

src/finiteArea/fields/areaFields/areaFields.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,8 @@ License
 
 #include "faMesh.H"
 #include "areaFields.H"
+#include "coupledFaPatchField.H"
+#include "registerSwitch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -46,6 +48,158 @@ defineTemplateTypeNameAndDebug(areaSphericalTensorField, 0);
 defineTemplateTypeNameAndDebug(areaSymmTensorField, 0);
 defineTemplateTypeNameAndDebug(areaTensorField, 0);
 
+defineTemplateDebugSwitchWithName
+(
+    areaScalarField::Boundary,
+    "areaScalarField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    areaVectorField::Boundary,
+    "areaVectorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    areaSphericalTensorField::Boundary,
+    "areaSphericalTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    areaSymmTensorField::Boundary,
+    "areaSymmTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    areaTensorField::Boundary,
+    "areaTensorField::Boundary",
+    0
+);
+
+template<> scalar areaScalarField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("areaScalarField::Boundary::tolerance", 0)
+);
+registerOptSwitch
+(
+    "areaScalarField::Boundary::tolerance",
+    scalar,
+    Foam::areaScalarField::Boundary::tolerance
+);
+
+template<> scalar areaVectorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("areaVectorField::Boundary::tolerance", 0)
+);
+registerOptSwitch
+(
+    "areaVectorField::Boundary::tolerance",
+    scalar,
+    Foam::areaVectorField::Boundary::tolerance
+);
+
+template<> scalar areaSphericalTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch
+    (
+        "areaSphericalTensorField::Boundary::tolerance",
+        0
+    )
+);
+registerOptSwitch
+(
+    "areaSphericalTensorField::Boundary::tolerance",
+    scalar,
+    Foam::areaSphericalTensorField::Boundary::tolerance
+);
+
+template<> scalar areaSymmTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch
+    (
+        "areaSymmTensorField::Boundary::tolerance",
+        0
+    )
+);
+registerOptSwitch
+(
+    "areaSymmTensorField::Boundary::tolerance",
+    scalar,
+    Foam::areaSymmTensorField::Boundary::tolerance
+);
+
+template<> scalar areaTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("areaTensorField::Boundary::tolerance", 0)
+);
+registerOptSwitch
+(
+    "areaTensorField::Boundary::tolerance",
+    scalar,
+    Foam::areaTensorField::Boundary::tolerance
+);
+
+
+// Local-ops consistency enforcing
+
+template<> int areaScalarField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "areaScalarField::Boundary::localConsistency",
+    int,
+    Foam::areaScalarField::Boundary::localConsistency
+);
+
+template<> int areaVectorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "areaVectorField::Boundary::localConsistency",
+    int,
+    Foam::areaVectorField::Boundary::localConsistency
+);
+
+template<> int areaSphericalTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "areaSphericalTensorField::Boundary::localConsistency",
+    int,
+    Foam::areaSphericalTensorField::Boundary::localConsistency
+);
+
+template<> int areaSymmTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "areaSymmTensorField::Boundary::localConsistency",
+    int,
+    Foam::areaSymmTensorField::Boundary::localConsistency
+);
+
+template<> int areaTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "areaTensorField::Boundary::localConsistency",
+    int,
+    Foam::areaTensorField::Boundary::localConsistency
+);
+
 } // End namespace Foam
 
 
@@ -76,6 +230,57 @@ void GeometricField<scalar, faPatchField, areaMesh>::replace
     *this == gsf;
 }
 
+template<>
+bool GeometricBoundaryField<scalar, faPatchField, areaMesh>::check() const
+{
+    return checkConsistency<coupledFaPatchField<scalar>>
+    (
+        Foam::areaScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
+template<>
+bool GeometricBoundaryField<vector, faPatchField, areaMesh>::check() const
+{
+    return checkConsistency<coupledFaPatchField<vector>>
+    (
+        Foam::areaScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
+template<>
+bool GeometricBoundaryField<sphericalTensor, faPatchField, areaMesh>::check
+() const
+{
+    return checkConsistency<coupledFaPatchField<sphericalTensor>>
+    (
+        Foam::areaScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
+template<>
+bool GeometricBoundaryField<symmTensor, faPatchField, areaMesh>::check() const
+{
+    return checkConsistency<coupledFaPatchField<symmTensor>>
+    (
+        Foam::areaScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
+template<>
+bool GeometricBoundaryField<tensor, faPatchField, areaMesh>::check() const
+{
+    return checkConsistency<coupledFaPatchField<tensor>>
+    (
+        Foam::areaScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
 } // End namespace Foam
 
 

src/finiteArea/fields/areaFields/areaFields.H

@@ -70,6 +70,17 @@ void GeometricField<scalar, faPatchField, areaMesh>::replace
     const GeometricField<scalar, faPatchField, areaMesh>& sf
 );
 
+template<>
+bool GeometricBoundaryField<scalar, faPatchField, areaMesh>::check() const;
+template<>
+bool GeometricBoundaryField<vector, faPatchField, areaMesh>::check() const;
+template<>
+bool GeometricBoundaryField<sphericalTensor, faPatchField, areaMesh>::check
+() const;
+template<>
+bool GeometricBoundaryField<symmTensor, faPatchField, areaMesh>::check() const;
+template<>
+bool GeometricBoundaryField<tensor, faPatchField, areaMesh>::check() const;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteArea/fields/edgeFields/edgeFields.C

@@ -27,6 +27,7 @@ License
 
 #include "faMesh.H"
 #include "edgeFields.H"
+#include "registerSwitch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -45,6 +46,115 @@ defineTemplateTypeNameAndDebug(edgeSphericalTensorField, 0);
 defineTemplateTypeNameAndDebug(edgeSymmTensorField, 0);
 defineTemplateTypeNameAndDebug(edgeTensorField, 0);
 
+defineTemplateDebugSwitchWithName
+(
+    edgeScalarField::Boundary,
+    "edgeScalarField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    edgeVectorField::Boundary,
+    "edgeVectorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    edgeSphericalTensorField::Boundary,
+    "edgeSphericalTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    edgeSymmTensorField::Boundary,
+    "edgeSymmTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    edgeTensorField::Boundary,
+    "edgeTensorField::Boundary",
+    0
+);
+
+template<> scalar edgeScalarField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("tolerance", 0)
+);
+template<> scalar edgeVectorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("tolerance", 0)
+);
+template<> scalar edgeSphericalTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("tolerance", 0)
+);
+template<> scalar edgeSymmTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("tolerance", 0)
+);
+template<> scalar edgeTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("tolerance", 0)
+);
+
+// Local-ops consistency enforcing
+
+template<> int edgeScalarField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "edgeScalarField::Boundary::localConsistency",
+    int,
+    Foam::edgeScalarField::Boundary::localConsistency
+);
+
+template<> int edgeVectorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "edgeVectorField::Boundary::localConsistency",
+    int,
+    Foam::edgeVectorField::Boundary::localConsistency
+);
+
+template<> int edgeSphericalTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "edgeSphericalTensorField::Boundary::localConsistency",
+    int,
+    Foam::edgeSphericalTensorField::Boundary::localConsistency
+);
+
+template<> int edgeSymmTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "edgeSymmTensorField::Boundary::localConsistency",
+    int,
+    Foam::edgeSymmTensorField::Boundary::localConsistency
+);
+
+template<> int edgeTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "edgeTensorField::Boundary::localConsistency",
+    int,
+    Foam::edgeTensorField::Boundary::localConsistency
+);
+
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteArea/fields/faPatchFields/basic/coupled/coupledFaPatchField.H

@@ -157,6 +157,27 @@ public:
                 const Pstream::commsTypes commsType
             );
 
+            //- Initialise the evaluation of the patch field after a local
+            //  operation
+            virtual void initEvaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {
+                initEvaluate(commsType);
+            }
+
+            //- Evaluate the patch field after a local operation (e.g. *=)
+            virtual void evaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {
+                evaluate(commsType);
+            }
+
             //- Return the matrix diagonal coefficients corresponding to the
             //  evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueInternalCoeffs

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.H

@@ -192,6 +192,24 @@ public:
         //- Evaluate the patch field
         virtual void evaluate(const Pstream::commsTypes commsType);
 
+        //- Initialise the evaluation of the patch field after a local
+        //  operation. Dummy since operating on a copy
+        virtual void initEvaluateLocal
+        (
+            const Pstream::commsTypes commsType =
+                Pstream::commsTypes::blocking
+        )
+        {}
+
+        //- Evaluate the patch field after a local operation (e.g. *=).
+        //  Dummy since operating on a copy
+        virtual void evaluateLocal
+        (
+            const Pstream::commsTypes commsType =
+                Pstream::commsTypes::blocking
+        )
+        {}
+
         //- Return patch-normal gradient
         virtual tmp<Field<Type>> snGrad() const;
 

src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.H

@@ -103,16 +103,6 @@ protected:
         //  Useful when initially constructed without a dictionary
         virtual void readDict(const dictionary& dict);
 
-        //- Set updated state
-        void setUpdated(bool state) noexcept
-        {
-            updated_ = state;
-        }
-
-        //- Set matrix manipulated state. Currently a no-op for faPatchField.
-        void setManipulated(bool state) noexcept
-        {}
-
 
 public:
 
@@ -224,11 +214,10 @@ public:
             return updated_;
         }
 
-        //- True if the matrix has already been manipulated.
+        //- Set updated state
-        //- Currently ignored (always false) for faPatchField
+        void setUpdated(bool state) noexcept
-        bool manipulatedMatrix() const noexcept
         {
-            return false;
+            updated_ = state;
         }
 
 
@@ -536,6 +525,23 @@ public:
                     Pstream::commsTypes::blocking
             );
 
+            //- Initialise the evaluation of the patch field after a local
+            //  operation
+            virtual void initEvaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
+
+            //- Evaluate the patch field after a local operation (e.g. *=)
+            virtual void evaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
+
             //- Return the matrix diagonal coefficients corresponding to the
             //- evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueInternalCoeffs

src/finiteArea/fields/faePatchFields/faePatchField/faePatchField.H

@@ -445,6 +445,26 @@ public:
         }
 
 
+        // Evaluation Functions
+
+            //- Initialise the evaluation of the patch field after a local
+            //  operation
+            virtual void initEvaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
+
+            //- Evaluate the patch field after a local operation (e.g. *=)
+            virtual void evaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
+
+
         // Mapping
 
         //- Map (and resize as needed) from self given a mapping object

src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C

@@ -151,7 +151,7 @@ Foam::SRF::SRFModel::Fcentrifugal() const
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            omega_ ^ (omega_ ^ (mesh_.C() - origin_))
+            omega_ ^ (omega_ ^ (mesh_.C().internalField() - origin_))
         )
     );
 }
@@ -182,7 +182,10 @@ Foam::vectorField Foam::SRF::SRFModel::velocity
 
 Foam::tmp<Foam::volVectorField> Foam::SRF::SRFModel::U() const
 {
-    return tmp<volVectorField>
+    const int oldLocal = volVectorField::Boundary::localConsistency;
+    volVectorField::Boundary::localConsistency = 0;
+    tmp<volVectorField> relPos(mesh_.C() - origin_);
+    tmp<volVectorField> tU
     (
         new volVectorField
         (
@@ -194,10 +197,12 @@ Foam::tmp<Foam::volVectorField> Foam::SRF::SRFModel::U() const
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            omega_
+            omega_ ^ (relPos - axis_*(axis_ & relPos))
-          ^ ((mesh_.C() - origin_) - axis_*(axis_ & (mesh_.C() - origin_)))
         )
     );
+    volVectorField::Boundary::localConsistency = oldLocal;
+
+    return tU;
 }
 
 

src/finiteVolume/cfdTools/general/include/gh.H

@@ -5,5 +5,10 @@
       ? g & (cmptMag(g.value())/mag(g.value()))*hRef
       : dimensionedScalar("ghRef", g.dimensions()*dimLength, 0)
     );
+    const int oldLocal = volVectorField::Boundary::localConsistency;
+    volVectorField::Boundary::localConsistency = 0;
+
     volScalarField gh("gh", (g & mesh.C()) - ghRef);
     surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+
+    volVectorField::Boundary::localConsistency = oldLocal;

src/finiteVolume/fields/fvPatchFields/basic/coupled/coupledFvPatchField.H

@@ -182,6 +182,27 @@ public:
                 const Pstream::commsTypes commsType
             );
 
+            //- Initialise the evaluation of the patch field after a local
+            //  operation
+            virtual void initEvaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {
+                initEvaluate(commsType);
+            }
+
+            //- Evaluate the patch field after a local operation (e.g. *=)
+            virtual void evaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {
+                evaluate(commsType);
+            }
+
             //- Return the matrix diagonal coefficients corresponding to the
             //  evaluation of the value of this patchField with given weights
             virtual tmp<Field<Type>> valueInternalCoeffs

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C

@@ -268,7 +268,13 @@ Foam::cyclicAMIFvPatchField<Type>::patchNeighbourField
     const Field<Type>& iField
 ) const
 {
-    // Bypass polyPatch to get nbrId. Instead use cyclicAMIFvPatch virtual
+    DebugPout
+        << "cyclicAMIFvPatchField::patchNeighbourField(const Field<Type>&) :"
+        << " field:" << this->internalField().name()
+        << " patch:" << this->patch().name()
+        << endl;
+
+    // By pass polyPatch to get nbrId. Instead use cyclicAMIFvPatch virtual
     // neighbPatch()
     const cyclicAMIFvPatch& neighbPatch = cyclicAMIPatch_.neighbPatch();
     const labelUList& nbrFaceCells = neighbPatch.faceCells();
@@ -500,6 +506,11 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
+    DebugPout<< "cyclicAMIFvPatchField::updateInterfaceMatrix() :"
+        << " field:" << this->internalField().name()
+        << " patch:" << this->patch().name()
+        << endl;
+
     const labelUList& faceCells = lduAddr.patchAddr(patchId);
 
     const auto& AMI =
@@ -618,6 +629,11 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
+    DebugPout<< "cyclicAMIFvPatchField::updateInterfaceMatrix() :"
+        << " field:" << this->internalField().name()
+        << " patch:" << this->patch().name()
+        << endl;
+
     const labelUList& faceCells = lduAddr.patchAddr(patchId);
 
     const auto& AMI = this->ownerAMI();

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H

@@ -200,6 +200,24 @@ public:
         //- Evaluate the patch field
         virtual void evaluate(const Pstream::commsTypes commsType);
 
+        //- Initialise the evaluation of the patch field after a local
+        //  operation. Dummy since operating on a copy
+        virtual void initEvaluateLocal
+        (
+            const Pstream::commsTypes commsType =
+                Pstream::commsTypes::blocking
+        )
+        {}
+
+        //- Evaluate the patch field after a local operation (e.g. *=).
+        //  Dummy since operating on a copy
+        virtual void evaluateLocal
+        (
+            const Pstream::commsTypes commsType =
+                Pstream::commsTypes::blocking
+        )
+        {}
+
         //- Return patch-normal gradient
         virtual tmp<Field<Type>> snGrad
         (

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H

@@ -102,7 +102,8 @@ class fvPatchFieldBase
         word patchType_;
 
 
-protected:
+//protected:
+public:
 
     // Protected Member Functions
 
@@ -613,6 +614,22 @@ public:
                     Pstream::commsTypes::blocking
             );
 
+            //- Initialise the evaluation of the patch field after a local
+            //  operation
+            virtual void initEvaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
+
+            //- Evaluate the patch field after a local operation (e.g. *=)
+            virtual void evaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
 
             //- Return the matrix diagonal coefficients corresponding to the
             //  evaluation of the value of this patchField with given weights

src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H

@@ -477,22 +477,39 @@ public:
 
         // Evaluation Functions
 
-        //- Initialise the evaluation of the patch field, generally a no-op
+            //- Initialise the evaluation of the patch field
             virtual void initEvaluate
             (
-            const Pstream::commsTypes commsType = Pstream::commsTypes::blocking
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
             )
             {}
 
-        //- Evaluate the patch field, generally a no-op
+            //- Evaluate the patch field, sets Updated to false
             virtual void evaluate
             (
-            const Pstream::commsTypes commsType = Pstream::commsTypes::blocking
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
             )
             {}
 
+            //- Initialise the evaluation of the patch field after a local
+            //  operation
+            virtual void initEvaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
+
+            //- Evaluate the patch field after a local operation (e.g. *=)
+            virtual void evaluateLocal
+            (
+                const Pstream::commsTypes commsType =
+                    Pstream::commsTypes::blocking
+            )
+            {}
 
-    // Other
 
         //- Write
         virtual void write(Ostream& os) const;

src/finiteVolume/fields/surfaceFields/surfaceFields.C

@@ -27,6 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "surfaceFields.H"
+#include "registerSwitch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -45,6 +46,125 @@ defineTemplateTypeNameAndDebug(surfaceSphericalTensorField, 0);
 defineTemplateTypeNameAndDebug(surfaceSymmTensorField, 0);
 defineTemplateTypeNameAndDebug(surfaceTensorField, 0);
 
+defineTemplateDebugSwitchWithName
+(
+    surfaceScalarField::Boundary,
+    "surfaceScalarField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    surfaceVectorField::Boundary,
+    "surfaceVectorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    surfaceSphericalTensorField::Boundary,
+    "surfaceSphericalTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    surfaceSymmTensorField::Boundary,
+    "surfaceSymmTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    surfaceTensorField::Boundary,
+    "surfaceTensorField::Boundary",
+    0
+);
+
+template<> scalar surfaceScalarField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("surfaceScalarField::Boundary::tolerance", 0)
+);
+template<> scalar surfaceVectorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("surfaceVectorField::Boundary::tolerance", 0)
+);
+template<> scalar surfaceSphericalTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch
+    (
+        "surfaceSphericalTensorField::Boundary::tolerance",
+        0
+    )
+);
+template<> scalar surfaceSymmTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch
+    (
+        "surfaceSymmTensorField::Boundary::tolerance",
+        0
+    )
+);
+template<> scalar surfaceTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("surfaceTensorField::Boundary::tolerance", 0)
+);
+
+
+// Local-ops consistency enforcing
+
+template<> int surfaceScalarField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "surfaceScalarField::Boundary::localConsistency",
+    int,
+    Foam::surfaceScalarField::Boundary::localConsistency
+);
+
+template<> int surfaceVectorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "surfaceVectorField::Boundary::localConsistency",
+    int,
+    Foam::surfaceVectorField::Boundary::localConsistency
+);
+
+template<> int surfaceSphericalTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "surfaceSphericalTensorField::Boundary::localConsistency",
+    int,
+    Foam::surfaceSphericalTensorField::Boundary::localConsistency
+);
+
+template<> int surfaceSymmTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "surfaceSymmTensorField::Boundary::localConsistency",
+    int,
+    Foam::surfaceSymmTensorField::Boundary::localConsistency
+);
+
+template<> int surfaceTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "surfaceTensorField::Boundary::localConsistency",
+    int,
+    Foam::surfaceTensorField::Boundary::localConsistency
+);
+
+
 } // End namespace Foam
 
 

src/finiteVolume/fields/volFields/volFields.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +27,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "volFields.H"
+#include "coupledFvPatchField.H"
+#include "registerSwitch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -45,6 +47,162 @@ defineTemplateTypeNameAndDebug(volSphericalTensorField, 0);
 defineTemplateTypeNameAndDebug(volSymmTensorField, 0);
 defineTemplateTypeNameAndDebug(volTensorField, 0);
 
+defineTemplateDebugSwitchWithName
+(
+    volScalarField::Boundary,
+    "volScalarField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    volVectorField::Boundary,
+    "volVectorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    volSphericalTensorField::Boundary,
+    "volSphericalTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    volSymmTensorField::Boundary,
+    "volSymmTensorField::Boundary",
+    0
+);
+defineTemplateDebugSwitchWithName
+(
+    volTensorField::Boundary,
+    "volTensorField::Boundary",
+    0
+);
+
+
+// Tolerance optimisation switch
+
+template<> scalar volScalarField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("volScalarField::Boundary::tolerance", 0)
+);
+registerOptSwitch
+(
+    "volScalarField::Boundary::tolerance",
+    scalar,
+    Foam::volScalarField::Boundary::tolerance
+);
+
+template<> scalar volVectorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("volScalarField::Boundary::tolerance", 0)
+);
+registerOptSwitch
+(
+    "volVectorField::Boundary::tolerance",
+    scalar,
+    Foam::volVectorField::Boundary::tolerance
+);
+
+template<> scalar volSphericalTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch
+    (
+        "volSphericalTensorField::Boundary::tolerance",
+        0
+    )
+);
+registerOptSwitch
+(
+    "volSphericalTensorField::Boundary::tolerance",
+    scalar,
+    Foam::volSphericalTensorField::Boundary::tolerance
+);
+
+template<> scalar volSymmTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch
+    (
+        "volSymmTensorField::Boundary::tolerance",
+        0
+    )
+);
+registerOptSwitch
+(
+    "volSymmTensorField::Boundary::tolerance",
+    scalar,
+    Foam::volSymmTensorField::Boundary::tolerance
+);
+
+template<> scalar volTensorField::Boundary::tolerance
+(
+    debug::floatOptimisationSwitch("volTensorField::Boundary::tolerance", 0)
+);
+registerOptSwitch
+(
+    "volTensorField::Boundary::tolerance",
+    scalar,
+    Foam::volTensorField::Boundary::tolerance
+);
+
+
+// Local-ops consistency enforcing
+
+template<> int volScalarField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "volScalarField::Boundary::localConsistency",
+    int,
+    Foam::volScalarField::Boundary::localConsistency
+);
+
+template<> int volVectorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "volVectorField::Boundary::localConsistency",
+    int,
+    Foam::volVectorField::Boundary::localConsistency
+);
+
+template<> int volSphericalTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "volSphericalTensorField::Boundary::localConsistency",
+    int,
+    Foam::volSphericalTensorField::Boundary::localConsistency
+);
+
+template<> int volSymmTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "volSymmTensorField::Boundary::localConsistency",
+    int,
+    Foam::volSymmTensorField::Boundary::localConsistency
+);
+
+template<> int volTensorField::Boundary::localConsistency
+(
+    debug::optimisationSwitch("localConsistency", 1)
+);
+registerOptSwitch
+(
+    "volTensorField::Boundary::localConsistency",
+    int,
+    Foam::volTensorField::Boundary::localConsistency
+);
+
+
 } // End namespace Foam
 
 
@@ -76,6 +234,60 @@ void GeometricField<scalar, fvPatchField, volMesh>::replace
     *this == gsf;
 }
 
+
+template<>
+bool GeometricBoundaryField<scalar, fvPatchField, volMesh>::check() const
+{
+    return checkConsistency<coupledFvPatchField<scalar>>
+    (
+        Foam::volScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
+template<>
+bool GeometricBoundaryField<vector, fvPatchField, volMesh>::check() const
+{
+    return checkConsistency<coupledFvPatchField<vector>>
+    (
+        Foam::volScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
+template<>
+bool GeometricBoundaryField<sphericalTensor, fvPatchField, volMesh>::check
+() const
+{
+    return checkConsistency<coupledFvPatchField<sphericalTensor>>
+    (
+        Foam::volScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
+template<>
+bool GeometricBoundaryField<symmTensor, fvPatchField, volMesh>::check
+() const
+{
+    return checkConsistency<coupledFvPatchField<symmTensor>>
+    (
+        Foam::volScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
+template<>
+bool GeometricBoundaryField<tensor, fvPatchField, volMesh>::check
+() const
+{
+    return checkConsistency<coupledFvPatchField<tensor>>
+    (
+        Foam::volScalarField::Boundary::tolerance,
+       !(debug&4)       // make into warning if debug&4
+    );
+}
+
 } // End namespace Foam
 
 

src/finiteVolume/fields/volFields/volFields.H

@@ -64,6 +64,18 @@ void GeometricField<scalar, fvPatchField, volMesh>::replace
     const GeometricField<scalar, fvPatchField, volMesh>& sf
 );
 
+template<>
+bool GeometricBoundaryField<scalar, fvPatchField, volMesh>::check() const;
+template<>
+bool GeometricBoundaryField<vector, fvPatchField, volMesh>::check() const;
+template<>
+bool GeometricBoundaryField<sphericalTensor, fvPatchField, volMesh>::check
+() const;
+template<>
+bool GeometricBoundaryField<symmTensor, fvPatchField, volMesh>::check() const;
+template<>
+bool GeometricBoundaryField<tensor, fvPatchField, volMesh>::check() const;
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/functionObjects/field/pressure/pressure.C

@@ -204,6 +204,9 @@ void Foam::functionObjects::pressure::addHydrostaticContribution
         (g_ & (cmptMag(g_.value())/mag(g_.value())))*hRef_
     );
 
+    const int oldLocal = volScalarField::Boundary::localConsistency;
+    volScalarField::Boundary::localConsistency = 0;
+
     tmp<volScalarField> rgh = rhoScale(p, (g_ & mesh_.C()) - ghRef);
 
     switch (hydrostaticMode_)
@@ -221,6 +224,8 @@ void Foam::functionObjects::pressure::addHydrostaticContribution
         default:
         {}
     }
+
+    volScalarField::Boundary::localConsistency = oldLocal;
 }
 
 

src/functionObjects/field/setFlow/setFlow.C

@@ -234,6 +234,9 @@ bool Foam::functionObjects::setFlow::execute()
         }
         case modeType::ROTATION:
         {
+            const int oldLocal = volVectorField::Boundary::localConsistency;
+            volVectorField::Boundary::localConsistency = 0;
+
             const volVectorField& C = mesh_.C();
             const volVectorField d
             (
@@ -266,6 +269,7 @@ bool Foam::functionObjects::setFlow::execute()
             }
 
             U = U & R_;
+            volVectorField::Boundary::localConsistency = oldLocal;
             U.correctBoundaryConditions();
             setPhi(U);
 
@@ -275,6 +279,9 @@ bool Foam::functionObjects::setFlow::execute()
         {
             const scalar pi = Foam::constant::mathematical::pi;
 
+            const int oldLocal = volVectorField::Boundary::localConsistency;
+            volVectorField::Boundary::localConsistency = 0;
+
             const volVectorField& C = mesh_.C();
 
             const volVectorField d
@@ -291,6 +298,7 @@ bool Foam::functionObjects::setFlow::execute()
             Uc.replace(vector::Z, sin(2*pi*x)*sqr(sin(pi*z)));
 
             U = U & R_;
+            volVectorField::Boundary::localConsistency = oldLocal;
             U.correctBoundaryConditions();
 
             // Calculating phi
@@ -334,6 +342,10 @@ bool Foam::functionObjects::setFlow::execute()
         case modeType::VORTEX3D:
         {
             const scalar pi = Foam::constant::mathematical::pi;
+
+            const int oldLocal = volVectorField::Boundary::localConsistency;
+            volVectorField::Boundary::localConsistency = 0;
+
             const volVectorField& C = mesh_.C();
 
             const volVectorField d
@@ -351,6 +363,7 @@ bool Foam::functionObjects::setFlow::execute()
             Uc.replace(vector::Z, -sin(2*pi*x)*sin(2*pi*y)*sqr(sin(pi*z)));
 
             U = U & R_;
+            volVectorField::Boundary::localConsistency = oldLocal;
             U.correctBoundaryConditions();
 
             // Calculating phi

src/functionObjects/field/surfaceDistance/surfaceDistance.C

@@ -144,7 +144,7 @@ bool Foam::functionObjects::surfaceDistance::execute()
 
     if (doCells_)
     {
-        const pointField& cc = mesh_.C();
+        const pointField& cc = mesh_.C().internalField();
 
         labelList surfaces;
         List<pointIndexHit> nearestInfo;
@@ -203,7 +203,7 @@ bool Foam::functionObjects::surfaceDistance::write()
 //
 //    if (doCells_)
 //    {
-//        const pointField& cc = mesh_.C();
+//        const pointField& cc = mesh_.C().internalField();
 //
 //        labelList surfaces;
 //        List<pointIndexHit> nearestInfo;