From adc38d55ac5e3eba776853d19e1f471277ac831e Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 24 May 2023 13:30:08 +0100
Subject: [PATCH] ENH: coupled: enforce consistency. See #2783

---
 .../turbulence/PDRkEpsilon/PDRkEpsilon.C      |  7 ++-
 .../liquidFilmFoam/liquidFilmFoam.C           |  8 ++-
 .../sphereSurfactantFoam/createFaFields.H     |  1 +
 .../preProcessing/setFields/setFields.C       | 13 ++++-
 .../fieldsDistributorTemplates.C              |  2 +-
 .../RAS/LamBremhorstKE/LamBremhorstKE.C       |  4 +-
 .../RAS/LienCubicKE/LienCubicKE.C             |  2 +
 .../RAS/LienLeschziner/LienLeschziner.C       |  4 +-
 .../RAS/ShihQuadraticKE/ShihQuadraticKE.C     |  2 +
 .../RAS/kkLOmega/kkLOmega.C                   |  2 +
 .../RAS/mixtureKEpsilon/mixtureKEpsilon.C     |  8 +++
 .../Base/kOmegaSST/kOmegaSSTBase.C            | 25 ++++++++++
 .../turbulenceModels/RAS/EBRSM/EBRSM.C        |  2 +
 .../turbulenceModels/RAS/LRR/LRR.C            |  2 +
 .../RAS/RNGkEpsilon/RNGkEpsilon.C             |  2 +
 .../turbulenceModels/RAS/SSG/SSG.C            |  2 +
 .../turbulenceModels/RAS/kEpsilon/kEpsilon.C  |  2 +
 .../RAS/kEpsilonPhitF/kEpsilonPhitF.C         |  4 +-
 .../turbulenceModels/RAS/kOmega/kOmega.C      |  2 +
 .../RAS/realizableKE/realizableKE.C           |  2 +
 .../ObukhovLength/ObukhovLength.C             | 50 ++++++++++++++++---
 .../kEpsilonLopesdaCosta.C                    |  2 +
 src/combustionModels/EDC/EDC.C                |  8 ++-
 src/combustionModels/PaSR/PaSR.C              |  5 +-
 .../cfdTools/general/bound/bound.C            |  4 ++
 .../CoEulerDdtScheme/CoEulerDdtScheme.C       | 42 ++++++++++++++--
 .../CrankNicolsonDdtScheme.C                  | 43 ++++++++++++++--
 .../EulerDdtScheme/EulerDdtScheme.C           | 37 ++++++++++++--
 .../ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C  | 42 ++++++++++++++--
 .../backwardDdtScheme/backwardDdtScheme.C     | 38 ++++++++++++--
 .../steadyStateDdtScheme.C                    | 28 +++++++----
 .../directionalMeshWavePatchDistMethod.C      |  6 ++-
 .../meshWave/meshWavePatchDistMethod.C        |  9 +++-
 .../meshWaveAddressingPatchDistMethod.C       | 22 ++++++++
 .../wallDistAddressing/wallDistAddressing.C   |  3 ++
 .../StandardChemistryModel.C                  |  3 +-
 36 files changed, 387 insertions(+), 51 deletions(-)

diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
index e8c36a390c0..c7152d8030f 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -132,6 +132,11 @@ void PDRkEpsilon::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push new cell values to
+    // coupled neighbours. Note that we want to avoid the re-updateCoeffs
+    // of the wallFunctions so make sure to bypass the evaluate on
+    // those patches and only do the coupled ones.
+    epsilon_.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
 
     // Add the blockage generation term so that it is included consistently
     // in both the k and epsilon equations
diff --git a/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C b/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C
index 4cf01b79a00..02878bf9e4f 100644
--- a/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C
+++ b/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C
@@ -87,7 +87,13 @@ int main(int argc, char *argv[])
             (
                 fam::ddt(h, Us)
               + fam::div(phi2s, Us)
-              + fam::Sp(0.0125*frictionFactor*mag(Us), Us)
+              + fam::Sp
+                (
+                    0.0125
+                   *frictionFactor.internalField()
+                   *mag(Us.internalField()),
+                    Us
+                )
              ==
                 Gs*h
               - fam::Sp(Sd, Us)
diff --git a/applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H b/applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H
index 8ba46024455..514c6157d40 100644
--- a/applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H
+++ b/applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H
@@ -47,6 +47,7 @@ forAll(Us, faceI)
     Us[faceI].z() =
         Uinf.value()*0.25*R[faceI].x()*R[faceI].z()/sqr(mag(R[faceI]));
 }
+Us.boundaryFieldRef().evaluateCoupled<coupledFaPatch>();
 
 Us -= aMesh.faceAreaNormals()*(aMesh.faceAreaNormals() & Us);
 
diff --git a/applications/utilities/preProcessing/setFields/setFields.C b/applications/utilities/preProcessing/setFields/setFields.C
index 2865c0e1df1..d4312085a79 100644
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +44,8 @@ Description
 #include "faceSet.H"
 #include "volFields.H"
 #include "areaFields.H"
+#include "coupledFvPatch.H"
+#include "coupledFaPatch.H"
 
 using namespace Foam;
 
@@ -186,6 +188,9 @@ bool setCellFieldType
             pfld = pfld.patchInternalField();
         }
 
+        // Handle any e.g. halo-swaps
+        field.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
+
         if (!field.write())
         {
             FatalErrorInFunction
@@ -271,6 +276,9 @@ bool setAreaFieldType
             }
         }
 
+        // Handle any e.g. halo-swaps
+        field.boundaryFieldRef().template evaluateCoupled<coupledFaPatch>();
+
         if (!field.write())
         {
             FatalErrorInFunction
@@ -421,6 +429,9 @@ bool setFaceFieldType
             }
         }
 
+        // Handle any e.g. halo-swaps
+        field.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
+
         if (!field.write())
         {
             FatalErrorInFunction
diff --git a/src/OpenFOAM/parallel/fieldsDistributor/fieldsDistributorTemplates.C b/src/OpenFOAM/parallel/fieldsDistributor/fieldsDistributorTemplates.C
index 31635263889..e260e16b53d 100644
--- a/src/OpenFOAM/parallel/fieldsDistributor/fieldsDistributorTemplates.C
+++ b/src/OpenFOAM/parallel/fieldsDistributor/fieldsDistributorTemplates.C
@@ -141,7 +141,7 @@ void Foam::fieldsDistributor::readFieldsImpl
         // only need to know about themselves.
         // Can broadcast decompose = yes/no from master
 
-        bitSet localValues(haveMeshOnProc)
+        bitSet localValues(haveMeshOnProc);
         bitSet masterValues(localValues);
         Pstream::broadcast(masterValues);
 
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LamBremhorstKE/LamBremhorstKE.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LamBremhorstKE/LamBremhorstKE.C
index 6844e216a47..30717fd6a3e 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LamBremhorstKE/LamBremhorstKE.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -221,6 +221,8 @@ void LamBremhorstKE::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
 
     const volScalarField Rt(this->Rt());
     const volScalarField fMu(this->fMu(Rt));
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKE/LienCubicKE.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKE/LienCubicKE.C
index cc03139c53a..9e25620447e 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKE/LienCubicKE.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKE/LienCubicKE.C
@@ -390,6 +390,8 @@ void LienCubicKE::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
 
     const volScalarField f2(this->f2());
 
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienLeschziner/LienLeschziner.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienLeschziner/LienLeschziner.C
index 21d0870960e..30a5c54dde8 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienLeschziner/LienLeschziner.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienLeschziner/LienLeschziner.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -274,6 +274,8 @@ void LienLeschziner::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
 
     const volScalarField f2(this->f2());
 
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
index 727ee0ddef8..1444ef3bb5c 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
@@ -274,6 +274,8 @@ void ShihQuadraticKE::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C
index 149940e33cf..e95b9b3ca4f 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C
@@ -679,6 +679,8 @@ void kkLOmega::correct()
 
 
     omega_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    omega_.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
 
     // Turbulence specific dissipation rate equation
     tmp<fvScalarMatrix> omegaEqn
diff --git a/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C b/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
index b2c743253ce..35aa10c5b07 100644
--- a/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
+++ b/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
@@ -619,7 +619,11 @@ void mixtureKEpsilon<BasicTurbulenceModel>::correct()
 
         // Update k, epsilon and G at the wall
         kl.boundaryFieldRef().updateCoeffs();
+        // Push any changed cell values to coupled neighbours
+        kl.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
+
         epsilonl.boundaryFieldRef().updateCoeffs();
+        epsilonl.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
 
         Gc.ref().checkOut();
     }
@@ -639,7 +643,11 @@ void mixtureKEpsilon<BasicTurbulenceModel>::correct()
 
         // Update k, epsilon and G at the wall
         kg.boundaryFieldRef().updateCoeffs();
+        // Push any changed cell values to coupled neighbours
+        kg.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
         epsilong.boundaryFieldRef().updateCoeffs();
+        // Push any changed cell values to coupled neighbours
+        epsilong.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
 
         Gd.ref().checkOut();
     }
diff --git a/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C
index 80f5bfa87af..ba6661b871d 100644
--- a/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C
+++ b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C
@@ -524,8 +524,33 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
     volScalarField::Internal GbyNu0(this->GbyNu0(tgradU(), S2));
     volScalarField::Internal G(this->GName(), nut*GbyNu0);
 
+
+    // - boundary condition changes a cell value
+    // - normally this would be triggered through correctBoundaryConditions
+    // - which would do
+    //      - fvPatchField::evaluate() which calls
+    //      - fvPatchField::updateCoeffs()
+    // - however any processor boundary conditions already start sending
+    //   at initEvaluate so would send over the old value.
+    // - avoid this by explicitly calling updateCoeffs early and then
+    //   only doing the boundary conditions that rely on initEvaluate
+    //   (currently only coupled ones)
+
+    //- 1. Explicitly swap values on coupled boundary conditions
     // Update omega and G at the wall
     omega_.boundaryFieldRef().updateCoeffs();
+    // omegaWallFunctions change the cell value! Make sure to push these to
+    // coupled neighbours. Note that we want to avoid the re-updateCoeffs
+    // of the wallFunctions so make sure to bypass the evaluate on
+    // those patches and only do the coupled ones.
+    omega_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
+
+    ////- 2. Make sure the boundary condition calls updateCoeffs from
+    ////     initEvaluate
+    ////     (so before any swap is done - requires all coupled bcs to be
+    ////      after wall bcs. Unfortunately this conflicts with cyclicACMI)
+    //omega_.correctBoundaryConditions();
+
 
     const volScalarField CDkOmega
     (
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/EBRSM/EBRSM.C b/src/TurbulenceModels/turbulenceModels/RAS/EBRSM/EBRSM.C
index db19a5deff0..cafc50905a6 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/EBRSM/EBRSM.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/EBRSM/EBRSM.C
@@ -474,6 +474,8 @@ void EBRSM<BasicTurbulenceModel>::correct()
 
         // Update epsilon and G at the wall
         epsilon_.boundaryFieldRef().updateCoeffs();
+        // Push any changed cell values to coupled neighbours
+        epsilon_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
         // (M:Eq. C.14)
         tmp<fvScalarMatrix> epsEqn
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/LRR/LRR.C b/src/TurbulenceModels/turbulenceModels/RAS/LRR/LRR.C
index cb20f989b88..9bd6a8d6660 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/LRR/LRR.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/LRR/LRR.C
@@ -295,6 +295,8 @@ void LRR<BasicTurbulenceModel>::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C
index c70797140b4..3f7df45d617 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -281,6 +281,8 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/SSG/SSG.C b/src/TurbulenceModels/turbulenceModels/RAS/SSG/SSG.C
index 721d7b86b58..6b1fa3291fe 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/SSG/SSG.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/SSG/SSG.C
@@ -305,6 +305,8 @@ void SSG<BasicTurbulenceModel>::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
index 62e3aa813ba..ead17ff2049 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
@@ -253,6 +253,8 @@ void kEpsilon<BasicTurbulenceModel>::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.C b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.C
index c3169d640ad..5e27438de2b 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilonPhitF/kEpsilonPhitF.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -412,6 +412,8 @@ void kEpsilonPhitF<BasicTurbulenceModel>::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
     // Turbulent kinetic energy dissipation rate equation (LUU:Eq. 4)
     // k/T ~ epsilon
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C
index 0035d5c666a..77721813cff 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C
@@ -211,6 +211,8 @@ void kOmega<BasicTurbulenceModel>::correct()
 
     // Update omega and G at the wall
     omega_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    omega_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
     // Turbulence specific dissipation rate equation
     tmp<fvScalarMatrix> omegaEqn
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.C b/src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.C
index 633ec4e00be..6e08625a61a 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.C
@@ -276,6 +276,8 @@ void realizableKE<BasicTurbulenceModel>::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
     // SAF: limiting thermo->nu(). If psiThermo is used rho might be < 0
     // temporarily when p < 0 then nu < 0 which needs limiting
diff --git a/src/atmosphericModels/functionObjects/ObukhovLength/ObukhovLength.C b/src/atmosphericModels/functionObjects/ObukhovLength/ObukhovLength.C
index 68e8d0eef7a..9c81e138140 100644
--- a/src/atmosphericModels/functionObjects/ObukhovLength/ObukhovLength.C
+++ b/src/atmosphericModels/functionObjects/ObukhovLength/ObukhovLength.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2020 ENERCON GmbH
-    Copyright (C) 2020 OpenCFD Ltd
+    Copyright (C) 2020,2023 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,6 +33,8 @@ License
 #include "mapPolyMesh.H"
 #include "addToRunTimeSelectionTable.H"
 #include "zeroGradientFvPatchField.H"
+#include "coupledFvPatchField.H"
+#include "processorFvPatchField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -125,16 +127,52 @@ bool Foam::functionObjects::ObukhovLength::calcOL()
             )
         );
 
+
+    // Special bit of coding here to handle cyclics
+    // - sign(B) can easily be positive in one cell but negative in its
+    //   coupled cell
+    // - so cyclic value might evaluate to 0
+    // - which upsets the division
+    // - note that instead of sign() any other field might have the same
+    //   positive and negative coupled values - it is just a lot less likely
+    // - problem is that overridden patch types do not propagate - use in-place
+    //   operations only
+
+    volScalarField denom
+    (
+        IOobject
+        (
+            "denom",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        ),
+        sign(B)
+    );
+
+    // Override interpolated value on interpolated coupled patches
+    for (auto& pfld : denom.boundaryFieldRef())
+    {
+        if
+        (
+            isA<coupledFvPatchField<scalar>>(pfld)
+        && !isA<processorFvPatchField<scalar>>(pfld)
+        )
+        {
+            pfld = pfld.patchInternalField();
+        }
+    }
+
+    denom *= kappa_*max(mag(B), dimensionedScalar(B.dimensions(), VSMALL));
+
     // (O:Eq. 26)
     *result1 = // ObukhovLength
         -min
         (
             dimensionedScalar(dimLength, ROOTVGREAT), // neutral stratification
-            pow3(*result2)/
-            (
-                sign(B)*kappa_
-                *max(mag(B), dimensionedScalar(B.dimensions(), VSMALL))
-            )
+            pow3(*result2)/denom
         );
 
     return isNew;
diff --git a/src/atmosphericModels/turbulenceModels/RAS/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C b/src/atmosphericModels/turbulenceModels/RAS/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C
index ffda9c451e1..7f6a4f8d676 100644
--- a/src/atmosphericModels/turbulenceModels/RAS/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C
+++ b/src/atmosphericModels/turbulenceModels/RAS/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C
@@ -421,6 +421,8 @@ void kEpsilonLopesdaCosta<BasicTurbulenceModel>::correct()
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
+    // Push any changed cell values to coupled neighbours
+    epsilon_.boundaryFieldRef().template evaluateCoupled<coupledFvPatch>();
 
     volScalarField::Internal magU(mag(U));
     volScalarField::Internal magU3(pow3(magU));
diff --git a/src/combustionModels/EDC/EDC.C b/src/combustionModels/EDC/EDC.C
index 5d22a9f137b..3e55ba9ce1c 100644
--- a/src/combustionModels/EDC/EDC.C
+++ b/src/combustionModels/EDC/EDC.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -143,6 +143,9 @@ void Foam::combustionModels::EDC<ReactionThermo>::correct()
                         );
                 }
             }
+
+            // Evaluate bcs
+            kappa_.correctBoundaryConditions();
         }
         else
         {
@@ -171,6 +174,9 @@ void Foam::combustionModels::EDC<ReactionThermo>::correct()
                         );
                 }
             }
+
+            // Evaluate bcs
+            kappa_.correctBoundaryConditions();
         }
 
         Info<< "Chemistry time solved max/min : "
diff --git a/src/combustionModels/PaSR/PaSR.C b/src/combustionModels/PaSR/PaSR.C
index 63a223c88d3..63559a0aafe 100644
--- a/src/combustionModels/PaSR/PaSR.C
+++ b/src/combustionModels/PaSR/PaSR.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -99,6 +99,9 @@ void Foam::combustionModels::PaSR<ReactionThermo>::correct()
                 kappa_[i] = 1.0;
             }
         }
+
+        // Evaluate bcs
+        kappa_.correctBoundaryConditions();
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/bound/bound.C b/src/finiteVolume/cfdTools/general/bound/bound.C
index 4e96f8c5a75..9c6c269ed56 100644
--- a/src/finiteVolume/cfdTools/general/bound/bound.C
+++ b/src/finiteVolume/cfdTools/general/bound/bound.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,6 +57,9 @@ Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
         );
 
         vsf.boundaryFieldRef() = max(vsf.boundaryField(), lowerBound.value());
+
+        // Give coupled bcs chance to update since cell values changed
+        vsf.boundaryFieldRef().evaluateCoupled<coupledFvPatch>();
     }
 
     return vsf;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index 2b16e0a36c5..14bd9c62178 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -178,6 +179,11 @@ CoEulerDdtScheme<Type>::fvcDdt
             rDeltaT.primitiveField()*dt.value()
            *(1.0 - mesh().Vsc0()/mesh().Vsc());
 
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
         return tdtdt;
     }
     else
@@ -214,7 +220,7 @@ CoEulerDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -232,6 +238,13 @@ CoEulerDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -266,7 +279,7 @@ CoEulerDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -284,6 +297,13 @@ CoEulerDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -318,7 +338,7 @@ CoEulerDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -339,6 +359,13 @@ CoEulerDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -374,7 +401,7 @@ CoEulerDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -404,6 +431,13 @@ CoEulerDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index a7ec7a31af4..3c395b5ae0f 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2021,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -393,6 +393,11 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
             (rDtCoef*dt)*(mesh().V() - mesh().V0())
           - mesh().V0()*offCentre_(ddt0.internalField())
         )/mesh().V();
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
     }
 
     return tdtdt;
@@ -447,7 +452,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
             );
         }
 
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -465,6 +470,13 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 ) - offCentre_(ff(ddt0.boundaryField()))
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -535,7 +547,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
             );
         }
 
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -555,6 +567,13 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 ) - offCentre_(ff(ddt0.boundaryField()))
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -629,7 +648,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
             );
         }
 
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -651,6 +670,13 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 ) - offCentre_(ff(ddt0.boundaryField()))
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -738,7 +764,7 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
             );
         }
 
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -772,6 +798,13 @@ CrankNicolsonDdtScheme<Type>::fvcDdt
                 ) - offCentre_(ff(ddt0.boundaryField()))
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index b94f628cad3..375c6daed8c 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -104,7 +105,7 @@ EulerDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -120,6 +121,13 @@ EulerDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -154,7 +162,7 @@ EulerDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -170,6 +178,13 @@ EulerDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -204,7 +219,7 @@ EulerDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -223,6 +238,13 @@ EulerDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -258,7 +280,7 @@ EulerDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -285,6 +307,13 @@ EulerDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
index 61e3186cdf3..9e2fcab6847 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -183,6 +184,11 @@ SLTSDdtScheme<Type>::fvcDdt
         tdtdt.ref().primitiveFieldRef() =
             rDeltaT.primitiveField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
 
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
         return tdtdt;
     }
     else
@@ -219,7 +225,7 @@ SLTSDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -237,6 +243,13 @@ SLTSDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -271,7 +284,7 @@ SLTSDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -289,6 +302,13 @@ SLTSDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -323,7 +343,7 @@ SLTSDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -344,6 +364,13 @@ SLTSDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -379,7 +406,7 @@ SLTSDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -409,6 +436,13 @@ SLTSDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index ef0f0aedf31..8d77b3c0095 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -152,7 +152,7 @@ backwardDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -178,6 +178,13 @@ backwardDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -224,7 +231,7 @@ backwardDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -250,6 +257,13 @@ backwardDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -296,7 +310,7 @@ backwardDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -325,6 +339,13 @@ backwardDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
@@ -372,7 +393,7 @@ backwardDdtScheme<Type>::fvcDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, fvPatchField, volMesh>>
+        tmp<GeometricField<Type, fvPatchField, volMesh>> tdtdt
         (
             new GeometricField<Type, fvPatchField, volMesh>
             (
@@ -418,6 +439,13 @@ backwardDdtScheme<Type>::fvcDdt
                 )
             )
         );
+
+        // Different operation on boundary v.s. internal so re-evaluate
+        // coupled boundaries
+        tdtdt.ref().boundaryFieldRef().
+            template evaluateCoupled<coupledFvPatch>();
+
+        return tdtdt;
     }
     else
     {
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
index a9702062487..849852f7544 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -374,20 +375,25 @@ tmp<surfaceScalarField> steadyStateDdtScheme<Type>::meshPhi
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
-    return tmp<surfaceScalarField>::New
+    auto tphi
     (
-        IOobject
+        tmp<surfaceScalarField>::New
         (
-            "meshPhi",
-            mesh().time().timeName(),
-            mesh().thisDb(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            IOobject::NO_REGISTER
-        ),
-        mesh(),
-        dimensionedScalar(dimVolume/dimTime, Zero)
+            IOobject
+            (
+                "meshPhi",
+                mesh().time().timeName(),
+                mesh().thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh(),
+            dimensionedScalar(dimVolume/dimTime, Zero)
+        )
     );
+    tphi.ref().setOriented();
+    return tphi;
 }
 
 
diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/directionalMeshWave/directionalMeshWavePatchDistMethod.C b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/directionalMeshWave/directionalMeshWavePatchDistMethod.C
index 45d36c44eed..7610085b70b 100644
--- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/directionalMeshWave/directionalMeshWavePatchDistMethod.C
+++ b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/directionalMeshWave/directionalMeshWavePatchDistMethod.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -150,6 +150,10 @@ bool Foam::patchDistMethods::directionalMeshWave::correct
         }
     }
 
+    // Make sure boundary values are up-to-date
+    y.correctBoundaryConditions();
+    n.correctBoundaryConditions();
+
     // Transfer number of unset values
     this->nUnset_ = wave.nUnset();
 
diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWave/meshWavePatchDistMethod.C b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWave/meshWavePatchDistMethod.C
index c3a13cf6338..0a4504b313d 100644
--- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWave/meshWavePatchDistMethod.C
+++ b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWave/meshWavePatchDistMethod.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -99,6 +99,9 @@ bool Foam::patchDistMethods::meshWave::correct(volScalarField& y)
         }
     }
 
+    // Make sure boundary values are up-to-date
+    y.correctBoundaryConditions();
+
     // Transfer number of unset values
     nUnset_ = wave.nUnset();
 
@@ -155,6 +158,10 @@ bool Foam::patchDistMethods::meshWave::correct
         }
     }
 
+    // Make sure boundary values are up-to-date
+    y.correctBoundaryConditions();
+    n.correctBoundaryConditions();
+
     // Transfer number of unset values
     nUnset_ = wave.nUnset();
 
diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWaveAddressing/meshWaveAddressingPatchDistMethod.C b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWaveAddressing/meshWaveAddressingPatchDistMethod.C
index bf9fb619b63..035173a4172 100644
--- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWaveAddressing/meshWaveAddressingPatchDistMethod.C
+++ b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/meshWaveAddressing/meshWaveAddressingPatchDistMethod.C
@@ -99,6 +99,26 @@ bool Foam::patchDistMethods::meshWaveAddressing::correct
 
     y = wDist.y();
 
+    // Note: copying value only so might not be consistent with supplied
+    // patch types (e.g. zeroGradient when called from wallDist). Assume
+    // only affected ones are the supplied patches ...
+    {
+        auto& bfld = y.boundaryFieldRef();
+        const label startOfRequests = UPstream::nRequests();
+        for (const label patchi : patchIDs_)
+        {
+            bfld[patchi].initEvaluate(UPstream::commsTypes::nonBlocking);
+        }
+
+        // Wait for outstanding requests
+        UPstream::waitRequests(startOfRequests);
+
+        for (const label patchi : patchIDs_)
+        {
+            bfld[patchi].evaluate(UPstream::commsTypes::nonBlocking);
+        }
+    }
+
 
     // Only calculate n if the field is defined
     if (notNull(n))
@@ -109,6 +129,8 @@ bool Foam::patchDistMethods::meshWaveAddressing::correct
             pnf == pnf.patch().nf();
         }
 
+        // No problem with inconsistency as for y (see above) since doing
+        // correctBoundaryConditions on actual n field.
         wDist.map(n, mapDistribute::transform());
     }
 
diff --git a/src/finiteVolume/fvMesh/wallDist/wallDistAddressing/wallDistAddressing.C b/src/finiteVolume/fvMesh/wallDist/wallDistAddressing/wallDistAddressing.C
index 80c97e1541c..c6639c28571 100644
--- a/src/finiteVolume/fvMesh/wallDist/wallDistAddressing/wallDistAddressing.C
+++ b/src/finiteVolume/fvMesh/wallDist/wallDistAddressing/wallDistAddressing.C
@@ -238,6 +238,9 @@ void Foam::wallDistAddressing::correct(volScalarField& y)
         );
     }
 
+    // Make sure boundary values are up to date
+    y.correctBoundaryConditions();
+
 
     // Extract all addressing
     // ~~~~~~~~~~~~~~~~~~~~~~
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
index 3b55b76d946..a306b2dadd4 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2021,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -574,6 +574,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::Qdot() const
                 Qdot[celli] -= hi*RR_[i][celli];
             }
         }
+        tQdot.ref().correctBoundaryConditions();
     }
 
     return tQdot;
-- 
GitLab