From 3e80552244ddbee03ad782a79a9a4b854b2be79c Mon Sep 17 00:00:00 2001
From: Johan Roenby <johan.roenby@gmail.com>
Date: Tue, 10 Dec 2024 14:29:56 +0000
Subject: [PATCH] BUG: isoAdvection: conserve volume of fluid across cyclic
patches (fixes #2457)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Co-authored-by: David Müller <>
Co-authored-by: Konstantinos Missios <>
---
.../isoAdvection/isoAdvection.C | 9 +++++
.../isoAdvection/isoAdvectionTemplates.C | 35 +++++++++++++++++++
2 files changed, 44 insertions(+)
diff --git a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
index a2296c3c6bf..6ecf2ad0176 100644
--- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
+++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
@@ -397,6 +397,15 @@ void Foam::isoAdvection::timeIntegratedFlux()
magSf
);
+ // Handling upwind cyclic boundary patches
+ const polyPatch& pp = boundaryMesh[patchi];
+ const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);
+ if (cpp)
+ {
+ const label neiPatchID(cpp->neighbPolyPatchID());
+ dVfb[neiPatchID][patchFacei] = -dVfb[patchi][patchFacei];
+ }
+
// Check if the face is on processor patch and append it to
// the list if necessary
checkIfOnProcPatch(facei);
diff --git a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
index c79592c3b4a..422d86c8963 100644
--- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
+++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
@@ -188,8 +188,43 @@ void Foam::isoAdvection::limitFluxes
+ faceValue(dVfcorrectionValues, facei);
setFaceValue(dVf_, facei, corrVf);
+
+ // If facei is on a cyclic patch correct dVf of facej on
+ // neighbour patch and alpha value of facej's owner cell.
+ if (!mesh_.isInternalFace(facei))
+ {
+ const polyBoundaryMesh& boundaryMesh =
+ mesh_.boundaryMesh();
+ const label patchi = boundaryMesh.patchID(facei);
+ const polyPatch& pp = boundaryMesh[patchi];
+ const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);
+ const label patchFacei = pp.whichFace(facei);
+
+ if (cpp)
+ {
+ const label neiPatchID(cpp->neighbPolyPatchID());
+ surfaceScalarField::Boundary& dVfb =
+ dVf_.boundaryFieldRef();
+ dVfb[neiPatchID][patchFacei] =
+ -dVfb[patchi][patchFacei];
+ const polyPatch& np = boundaryMesh[neiPatchID];
+ const label globalFacei = np.start() + patchFacei;
+ const label neiOwn(owner[globalFacei]);
+ scalar VneiOwn = mesh_.V()[neiOwn];
+ if (porosityEnabled_)
+ {
+ VneiOwn *=
+ porosityPtr_->primitiveField()[neiOwn];
+ }
+ alpha1_[neiOwn] +=
+ faceValue(dVfcorrectionValues, facei)/VneiOwn;
+ }
+ }
+
}
+
}
+
syncProcPatches(dVf_, phi_);
}
else
--
GitLab