From 4fe3f55e4d2cf7250ae17653fc26ba2f0f0666ad Mon Sep 17 00:00:00 2001
From: Johan Roenby <johan.roenby@gmail.com>
Date: Fri, 10 Jan 2025 15:32:29 +0000
Subject: [PATCH] BUG: Fixes bug in plicRDF interface reconstruction across
 cyclic boundaries (fixes #3279)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Description:
The plicRDF interface reconstruction scheme calculates a reconstructed distance
function (RDF) in all interface cells and their point neighbours.
In point neighbours to interface cells, the RDF is calculated as a weighted
average of the distances to all adjacent interface cell centres with the weight
being the inverse distance to the cell centre.
By using the zoneDistribute class written by Henning Scheufler, the required
stencil data is communicated efficiently for stencil cells living on different
sides of one or more processor patches.
Some of the data required for the RDF reconstruction are cell centre and
interface centre positions. When a stencil extends across a cyclic patch these
positions have so far not been properly transformed in OpenFOAM. This issue is
fixed by the current contribution.
The fix is done by modifying the zoneDistribute class to hold the required
information about zone cells adjacent to cyclic patches. Positions are then
communicated with a new getPosition function which replaces getValue for
position data in the reconstructedDistanceFunction and plicRDF classes.
The implementation does not change the behaviour for cells not on
a cyclic patch and should have insignificant effect on efficiency for these.
The implementation can probably be optimised in terms of efficiency for zone
cells on cyclic patches, but we note that there are generally only very few of
these (interface cells and their point neighbours on cyclic patches) and so
the potential for speedup is expected to be limited.

Current limitations:

- In parallel, the user must constrain the decomposition to preserve cyclic
patches on the same processor, for the implementation to work properly.
  - See an example here: tutorials/discInConstantFlowCyclicBCs/system/decomposeParDict
- In the case of parallel rotational cyclics that are split by the decomposition
the current bugfix does not work and therefore throws an error. This is ongoing
work and should be reported and fixed by a future patch.

For further details, please see the modified files and the comments therein:
- $FOAM_SRC/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
- $FOAM_SRC/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
- $FOAM_SRC/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
- $FOAM_SRC/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
- $FOAM_SRC/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C

Co-authored-by: David Müller <> KIT
Co-authored-by: Konstantinos Missios <> Roskilde University
Co-authored-by: Johan Roenby <> Roskilde University and STROMNING
---
 .../fvMesh/zoneDistribute/zoneDistribute.C    | 262 +++++++++++++++++-
 .../fvMesh/zoneDistribute/zoneDistribute.H    |  41 +++
 .../fvMesh/zoneDistribute/zoneDistributeI.H   |  54 +++-
 .../reconstructedDistanceFunction.C           |  31 ++-
 .../plicSchemes/plicRDF/plicRDF.C             |  24 +-
 .../discInConstantFlowCyclicBCs/Allclean      |   2 +
 .../discInConstantFlowCyclicBCs/Allrun        |  12 +-
 .../system/controlDict                        |   2 +-
 .../system/decomposeParDict                   |   8 +
 .../system/fvSolution                         |   6 +-
 10 files changed, 424 insertions(+), 18 deletions(-)

diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
index 0aa301237d..b52943570c 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
@@ -27,6 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "zoneDistribute.H"
+#include "processorPolyPatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -44,10 +45,48 @@ Foam::zoneDistribute::zoneDistribute(const fvMesh& mesh)
     stencil_(zoneCPCStencil::New(mesh)),
     globalNumbering_(stencil_.globalNumbering()),
     send_(UPstream::nProcs()),
-    pBufs_(UPstream::commsTypes::nonBlocking)
+    pBufs_(UPstream::commsTypes::nonBlocking),
+    cyclicBoundaryCells_(mesh.nCells(), false)
 {
     // Don't clear storage on persistent buffer
     pBufs_.allowClearRecv(false);
+
+    // Loop over boundary patches and store cells with a face on a cyclic patch
+    bool hasCyclicPatches = false;
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const cyclicPolyPatch* cpp =
+            isA<cyclicPolyPatch>(mesh.boundaryMesh()[patchi]);
+
+        if (cpp)
+        {
+            cyclicBoundaryCells_.set(cpp->faceCells());
+            hasCyclicPatches = true;
+        }
+    }
+
+    // Populate cyclicCentres_
+    if(hasCyclicPatches)
+    {
+        // Make a boolList from the bitSet
+        boolList isCyclicCell(mesh.nCells(), false);
+
+        forAll(cyclicBoundaryCells_, celli)
+        {
+            if (cyclicBoundaryCells_.test(celli))
+            {
+                isCyclicCell[celli] = true;
+            }
+        }
+
+        // Use getFields to get map of cell centres across processor boundaries
+        setUpCommforZone(isCyclicCell, true);
+
+        cyclicCentres_.reset
+        (
+            new Map<vectorField>(getFields(isCyclicCell, mesh_.C()))
+        );
+    }
 }
 
 
@@ -141,5 +180,226 @@ void Foam::zoneDistribute::setUpCommforZone
     }
 }
 
+Foam::List<Foam::label> Foam::zoneDistribute::getCyclicPatches
+(
+    const label celli,
+    const label globalIdx
+) const
+{
+    // Initialise cyclic patch label list
+    List<label> patches(0);
+
+    // If celli is not on a cyclic boundary, return the empty list
+    if (!cyclicBoundaryCells_.test(celli))
+    {
+        return patches;
+    }
+
+    const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
+
+    // So celli belongs to at least one cyclic patch.
+    // Let us figure out which.
+    if (globalNumbering_.isLocal(globalIdx)) // celli and globalIdx on same proc
+    {
+        // Making list of cyclic patches to which celli belongs
+        List<label> celliCyclicPatches;
+        forAll(bMesh, patchi)
+        {
+            if (isA<cyclicPolyPatch>(bMesh[patchi]))
+            {
+                // Note: Probably not efficient due to use of found(celli) but
+                // typically only used for very cells (interface cells and their
+                // point neighbours on cyclic boundaries).
+                if (bMesh[patchi].faceCells().found(celli))
+                {
+                    celliCyclicPatches.append(patchi);
+                }
+            }
+        }
+
+        // Get all local point neighbor cells of celli, i.e. all point
+        // neighbours that are not on the other side of a cyclic patch.
+        List<label> localPointNeiCells(0);
+        const labelList& cellPoints = mesh_.cellPoints()[celli];
+
+        for (const label cellPoint : cellPoints)
+        {
+            const labelList& pointKCells = mesh_.pointCells()[cellPoint];
+
+            for (const label pointKCell : pointKCells)
+            {
+                if (!localPointNeiCells.found(pointKCell))
+                {
+                    localPointNeiCells.append(pointKCell);
+                }
+            }
+        }
+
+        // Since globalIdx is a global cell index obtained from the point
+        // neighbour list, stencil[celli], all cells in this that are not in
+        // localPointNeiCells must be cyclic neighbour cells.
+        const label localIdx = globalNumbering_.toLocal(globalIdx);
+        if (!localPointNeiCells.found(localIdx))
+        {
+            for (const label patchi : celliCyclicPatches)
+            {
+                // find the corresponding cyclic neighbor patch ID
+                const cyclicPolyPatch& cpp =
+                    static_cast<const cyclicPolyPatch&>(bMesh[patchi]);
+
+                const label neiPatch = cpp.neighbPatchID();
+
+                // Check if the cell globalIdx is on neiPatch.
+                // If it is, append neiPatch to list of
+                if (bMesh[neiPatch].faceCells().found(localIdx))
+                {
+                    patches.append(neiPatch);
+                    // Here it may be possible to append patchi and do:
+                    //
+                    //    cpp.transformPosition()
+                    //
+                    // instead of
+                    //
+                    //    cpp.neighbPatch().transformPosition() in getPosition()
+                }
+            }
+        }
+    }
+    else // celli and globalIdx on differet processors
+    {
+        // Note: The following is needed if a celli is located at the interface
+        // (plicRDF), on a cyclic patch and a processor patch. In this case
+        // globalIdx may be on a different processor, but requires the
+        // transformation from a cyclic patch on the processor of celli.
+
+        const List<label>& faces = mesh_.cells()[celli];
+
+        // Loop over all faces of celli and find cyclic patches
+        for (const label facei : faces)
+        {
+            if (mesh_.isInternalFace(facei)) continue;
+
+            const label patchi = bMesh.whichPatch(facei);
+            const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(bMesh[patchi]);
+
+            if (cpp)
+            {
+                // Get the neighbor cell across the cyclic face
+                const label cycNeiPatch = cpp->neighbPatchID();
+                const label cycNeiCell =
+                    bMesh[cycNeiPatch].faceCells()[facei - cpp->start()];
+                const List<label>& cycNeiCellFaces = mesh_.cells()[cycNeiCell];
+
+                // Loop over all the faces of the neighbor cell
+                for (const label cycNeiCellFace : cycNeiCellFaces)
+                {
+                    if (mesh_.isInternalFace(cycNeiCellFace))
+                    {
+                        continue;
+                    }
+
+                    // Check if the neighbor cell has processor patches
+                    const label neiPatch = bMesh.whichPatch(cycNeiCellFace);
+                    const processorPolyPatch* ppp =
+                        isA<processorPolyPatch>(bMesh[neiPatch]);
+
+                    if (ppp)
+                    {
+                        // Avoid duplicate entries
+                        if (patches.found(cycNeiPatch))
+                        {
+                            continue;
+                        }
+
+                        // Since we can not access any information on globalIdx
+                        // we use the cell centre map from the stencil to
+                        // identify the cell.
+                        const label localFaceID = cycNeiCellFace - ppp->start();
+                        const vector neiCentre =
+                            ppp->neighbFaceCellCentres()[localFaceID];
+
+                        forAll(cyclicCentres_()[celli], k)
+                        {
+                            if
+                            (
+                                (
+                                    mag(cyclicCentres_()[celli][k] - neiCentre)
+                                        < 100*SMALL
+                                )
+                                && (stencil_[celli][k] == globalIdx)
+                            )
+                            {
+                                patches.append(cycNeiPatch);
+                                // Here an alternative might be to append patchi
+                                // and do:
+                                //     cpp.transformPosition()
+                                // instead of
+                                //     cpp.neighbPatch().transformPosition()
+                                // in getValue()
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    return patches;
+}
+
+
+Foam::vector Foam::zoneDistribute::getPosition
+(
+    const VolumeField<vector>& positions,
+    const Map<vector>& valuesFromOtherProc,
+    const label gblIdx,
+    const List<label> cyclicPatchID
+) const
+{
+    // Position vector, possibly from other processor, to be returned
+    vector position(getValue(positions, valuesFromOtherProc, gblIdx));
+
+    // Dealing with position transformation across cyclic patches.
+    // If no transformation is required (most cases), cyclicPatchID is empty
+    forAll(cyclicPatchID, i)
+    {
+        const label patchi = cyclicPatchID[i];
+
+        const cyclicPolyPatch& cpp =
+            static_cast<const cyclicPolyPatch&>
+                (
+                    positions.mesh().boundaryMesh()[patchi]
+                );
+
+        if (cpp.transform() != coupledPolyPatch::transformType::ROTATIONAL)
+        {
+            cpp.neighbPatch().transformPosition(position, 0);
+        }
+        else if (globalNumbering_.isLocal(gblIdx))
+        {
+            const label localIdx = globalNumbering_.toLocal(gblIdx);
+
+            for (const label facei : mesh_.cells()[localIdx])
+            {
+                if (mesh_.boundaryMesh().whichPatch(facei) == cyclicPatchID[i])
+                {
+                    cpp.neighbPatch().transformPosition(position, facei);
+                    continue;
+                }
+            }
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Rotational cyclic patches are not supported in parallel.\n"
+                << "Try to decompose the domain so that the rotational cyclic"
+                << "patch is not split in between processors."
+                << exit(FatalError);
+        }
+    }
+
+    return position;
+}
+
 
 // ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
index dce29d9f86..575fcbd51e 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
@@ -115,6 +115,15 @@ class zoneDistribute
         //- Persistent set of exchange buffers
         PstreamBuffers pBufs_;
 
+        //- Cell labels for all cells with a face on a cyclic boundary
+        bitSet cyclicBoundaryCells_;
+
+        //- Holds for each cell on cyclic patch the centres of the cells in its
+        //  point neighbour stencil.
+        //  Used in getCyclicPatches to identify neigbhour patch ID of point
+        //  neighbours on other side of a processorPolyPatch.
+        autoPtr<Map<vectorField>> cyclicCentres_;
+
 
     // Private Member Functions
 
@@ -175,6 +184,19 @@ public:
             return globalNumbering_;
         }
 
+        //- Finds and returns list of all cyclic patch labels to which celli's
+        //  point neighbour cell, globalIdx, belongs. celli and globalIdx touch
+        //  in at least one point on these patces. globalIdx typically belongs
+        //  to stencil_[celli]. The returned label list is used to transform
+        //  positions across cyclic boundaries e.g. to be able to calculate
+        //  distances between cell centres and interface centres in plicRDF
+        //  across such boundaries.
+        List<label> getCyclicPatches
+        (
+            const label celli,
+            const label globalIdx
+        ) const;
+
         //- Gives patchNumber and patchFaceNumber for a given
         //- Geometric volume field
         template<typename Type>
@@ -185,6 +207,14 @@ public:
             const label gblIdx
         ) const;
 
+        vector getPosition
+        (
+            const VolumeField<vector>& positions,
+            const Map<vector>& valuesFromOtherProc,
+            const label gblIdx,
+            const List<label> cyclicPatchID = List<label>()
+        ) const;
+
         //- Returns stencil and provides a Map with globalNumbering
         template<typename Type>
         Map<Field<Type>> getFields
@@ -193,6 +223,17 @@ public:
             const VolumeField<Type>& phi
         );
 
+        //- Returns stencil and provides a Map with globalNumbering
+        //  Note: Not used currently (v2412) but needed for future surface
+        //  tension modelling.
+        template<typename Type>
+        Map<Field<Type>> getPositionFields
+        (
+            const boolList& zone,
+            const VolumeField<Type>& phi,
+            const bool& checkTransformation = false
+        );
+
         //- Returns stencil and provides a Map with globalNumbering
         template<typename Type>
         Map<Type> getDatafromOtherProc
diff --git a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
index 9be368d10f..e94beb1719 100644
--- a/src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
+++ b/src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
@@ -128,7 +128,59 @@ Foam::Map<Foam::Field<Type>> Foam::zoneDistribute::getFields
 
             for (const label gblIdx : stencil_[celli])
             {
-                tmpField.append(getValue(phi,neiValues,gblIdx));
+                tmpField.append(getValue(phi, neiValues, gblIdx));
+            }
+
+            stencilWithValues.emplace(celli, tmpField);
+        }
+    }
+
+    return stencilWithValues;
+}
+
+template<typename Type>
+Foam::Map<Foam::Field<Type>> Foam::zoneDistribute::getPositionFields
+(
+    const boolList& zone,
+    const VolumeField<Type>& phi,
+    const bool& checkTransformation
+)
+{
+    if (zone.size() != phi.size())
+    {
+        FatalErrorInFunction
+            << "size of zone: " << zone.size()
+            << "size of phi:" <<  phi.size()
+            << "do not match. Did the mesh change?"
+            << exit(FatalError);
+    }
+
+
+    // Get values from other proc
+    Map<Type> neiValues = getDatafromOtherProc(zone, phi);
+
+    Map<Field<Type>> stencilWithValues;
+
+    DynamicField<Type> tmpField(128);
+
+    forAll(zone, celli)
+    {
+        if (zone[celli])
+        {
+            tmpField.clear();
+
+            for (const label gblIdx : stencil_[celli])
+            {
+                List<label> cyclicPatches(0);
+                if(checkTransformation)
+                {
+                    cyclicPatches = getCyclicPatches(celli, gblIdx);
+                }
+
+                tmpField.append
+                (
+                    getPosition(phi, neiValues, gblIdx, cyclicPatches)
+                );
             }
 
             stencilWithValues.emplace(celli, tmpField);
diff --git a/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C b/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
index da0c0a0b47..5b0cd84e63 100644
--- a/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
+++ b/src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
@@ -282,9 +282,18 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
                     if (mag(n) != 0)
                     {
                         n /= mag(n);
-                        vector c = distribute.getValue(centre,mapCentres,gblIdx);
-                        vector distanceToIntSeg = (c - p);
-                        scalar distToSurf = distanceToIntSeg & (n);
+                        vector c
+                        (
+                            distribute.getPosition
+                            (
+                                centre,
+                                mapCentres,
+                                gblIdx,
+                                distribute.getCyclicPatches(celli, gblIdx)
+                            )
+                        );
+                        vector distanceToIntSeg(c - p);
+                        scalar distToSurf = distanceToIntSeg & n;
                         scalar weight = 0;
 
                         if (mag(distanceToIntSeg) != 0)
@@ -336,10 +345,18 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
                         if (mag(n) != 0)
                         {
                             n /= mag(n);
-                            vector c =
-                                distribute.getValue(centre, mapCentres, gblIdx);
-                            vector distanceToIntSeg = (c - p);
-                            scalar distToSurf = distanceToIntSeg & (n);
+                            vector c
+                            (
+                                distribute.getPosition
+                                (
+                                    centre,
+                                    mapCentres,
+                                    gblIdx,
+                                    distribute.getCyclicPatches(pCellI, gblIdx)
+                                )
+                            );
+                            vector distanceToIntSeg(c - p);
+                            scalar distToSurf = distanceToIntSeg & n;
                             scalar weight = 0;
 
                             if (mag(distanceToIntSeg) != 0)
diff --git a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C
index 62db3416f8..8d591f2371 100644
--- a/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C
+++ b/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C
@@ -96,8 +96,16 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
             if (mag(n) != 0)
             {
                 n /= mag(n);
-                vector centre =
-                    exchangeFields.getValue(centre_, mapCentre, gblIdx);
+                vector centre
+                (
+                    exchangeFields.getPosition
+                    (
+                        centre_,
+                        mapCentre,
+                        gblIdx,
+                        exchangeFields.getCyclicPatches(celli, gblIdx)
+                    )
+                );
                 vector distanceToIntSeg = (tensor::I- n*n) & (p - centre);
                 estimatedNormal += n /max(mag(distanceToIntSeg), SMALL);
                 weight += 1/max(mag(distanceToIntSeg), SMALL);
@@ -146,7 +154,11 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
                 const label gblIdx = stencil[celli][i];
                 cellCentre.append
                 (
-                    exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
+                    exchangeFields.getPosition
+                    (
+                        mesh_.C(), mapCC, gblIdx,
+                        exchangeFields.getCyclicPatches(celli, gblIdx)
+                    )
                 );
                 alphaValues.append
                 (
@@ -193,7 +205,11 @@ void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
         {
             cellCentre.append
             (
-                exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
+                exchangeFields.getPosition
+                (
+                    mesh_.C(), mapCC, gblIdx,
+                    exchangeFields.getCyclicPatches(celli, gblIdx)
+                )
             );
             phiValues.append
             (
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/Allclean b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/Allclean
index fb1f384730..c13797ec72 100755
--- a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/Allclean
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/Allclean
@@ -5,4 +5,6 @@ cd "${0%/*}" || exit                                # Run from this directory
 
 cleanCase0
 
+rm -f Phase1
+
 #------------------------------------------------------------------------------
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/Allrun b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/Allrun
index 321d8db8b1..8070263534 100755
--- a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/Allrun
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/Allrun
@@ -5,10 +5,20 @@ cd "${0%/*}" || exit                                # Run from this directory
 
 restore0Dir
 
+touch case.foam
+
 runApplication blockMesh
 
 runApplication setAlphaField
 
-runApplication $(getApplication)
+runApplication decomposePar
+
+runParallel $(getApplication)
+
+grep "Phase-1" log.$(getApplication) | cut -d' ' -f5 > Phase1
+
+echo "Volume conservation:"
+echo "Initial Phase-1: " $(head -1 Phase1)
+echo "Final Phase-1  : " $(tail -1 Phase1)
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/controlDict b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/controlDict
index bfd00f8afb..beccd14af8 100644
--- a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/controlDict
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/controlDict
@@ -35,7 +35,7 @@ purgeWrite      0;
 
 writeFormat     ascii;
 
-writePrecision  6;
+writePrecision  12;
 
 writeCompression off;
 
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/decomposeParDict b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/decomposeParDict
index 6c6f4b95f8..0caa1f4b9e 100644
--- a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/decomposeParDict
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/decomposeParDict
@@ -24,5 +24,13 @@ coeffs
     n   ( 2 1 2 );
 }
 
+constraints
+{
+    preservePatches
+    {
+        type            preservePatches;
+        patches         (left right top bottom);
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/fvSolution b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/fvSolution
index f242d11157..9a18a72b82 100644
--- a/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/fvSolution
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlowCyclicBCs/system/fvSolution
@@ -22,10 +22,10 @@ solvers
     {
         isoFaceTol      1e-6;
         surfCellTol     1e-6;
-        snapTol         1e-10;
+        snapTol         0;
         nAlphaBounds    3;
-        clip            true;
-        reconstructionScheme isoAlpha; // isoAlpha
+        clip            false;
+        reconstructionScheme plicRDF;
 
         nAlphaSubCycles 1;
         cAlpha          1; // Read but not used by interIsoFoam