INT: Integration of isoAdvector and supporting material

Community contribution from Johan Roenby, DHI

IsoAdvector is a geometric Volume-of-Fluid method for advection of a
sharp interface between two incompressible fluids. It works on both
structured and unstructured meshes with no requirements on cell shapes.
IsoAdvector is as an alternative choice for the interface compression
treatment with the MULES limiter implemented in the interFoam family
of solvers.

The isoAdvector concept and code was developed at DHI and was funded
by a Sapere Aude postdoc grant to Johan Roenby from The Danish Council
for Independent Research | Technology and Production Sciences (Grant-ID:
DFF - 1337-00118B - FTP).
Co-funding is also provided by the GTS grant to DHI from the Danish
Agency for Science, Technology and Innovation.

The ideas behind and performance of the isoAdvector scheme is
documented in:

    Roenby J, Bredmose H, Jasak H. 2016 A computational method for sharp
    interface  advection. R. Soc. open sci. 3: 160405.
    [http://dx.doi.org/10.1098/rsos.160405](http://dx.doi.org/10.1098/rsos.160405)

Videos showing isoAdvector's performance with a number of standard
test cases can be found in this youtube channel:

    https://www.youtube.com/channel/UCt6Idpv4C8TTgz1iUX0prAA

Project contributors:

* Johan Roenby <jro@dhigroup.com> (Inventor and main developer)
* Hrvoje Jasak <hrvoje.jasak@fsb.hr> (Consistent treatment of
  boundary faces including processor boundaries, parallelisation,
  code clean up
* Henrik Bredmose <hbre@dtu.dk> (Assisted in the conceptual
  development)
* Vuko Vukcevic <vuko.vukcevic@fsb.hr> (Code review, profiling,
  porting to foam-extend, bug fixing, testing)
* Tomislav Maric <tomislav@sourceflux.de> (Source file
  rearrangement)
* Andy Heather <a.heather@opencfd.co.uk> (Integration into OpenFOAM
  for v1706 release)

See the integration repository below to see the full set of changes
implemented for release into OpenFOAM v1706

    https://develop.openfoam.com/Community/Integration-isoAdvector

applications/solvers/multiphase/interIsoFoam/Make/files

+interIsoFoam.C
+
+EXE = $(FOAM_APPBIN)/interIsoFoam

applications/solvers/multiphase/interIsoFoam/Make/options

+EXE_INC = \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling

applications/solvers/multiphase/interIsoFoam/UEqn.H

+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(rho, U)
+     ==
+        fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+    }

applications/solvers/multiphase/interIsoFoam/alphaControls.H

+const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+
+label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+

applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    alphaCourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        mixture.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/interIsoFoam/alphaEqn.H

+// If there are more than one outer corrector, we use a mixture of old and
+// new U and phi for propagating alpha1 in all but the first outer iteration
+if (!pimple.firstIter())
+{
+    // We are recalculating alpha1 from the its old time value
+    alpha1 = alpha1.oldTime();
+    // Temporarily storing new U and phi values in prevIter storage
+    U.storePrevIter();
+    phi.storePrevIter();
+
+    // Overwriting new U and phi values with mixture of old and new values
+    phi = 0.5*(phi + phi.oldTime());
+    U = 0.5*(U + U.oldTime());
+}
+
+// Update alpha1
+advector.advect();
+
+// Update rhoPhi
+rhoPhi = advector.getRhoPhi(rho1, rho2);
+
+alpha2 = 1.0 - alpha1;
+
+if (!pimple.firstIter())
+{
+    // Restoring new U and phi values temporarily saved in prevIter() above
+    U = U.prevIter();
+    phi = phi.prevIter();
+}
+
+Info<< "Phase-1 volume fraction = "
+    << alpha1.weightedAverage(mesh.Vsc()).value()
+    << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+    << "  Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
+    << endl;

applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H

+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", rhoPhi.dimensions(), 0)
+    );
+
+    tmp<volScalarField> trSubDeltaT;
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + alpha2*rho2;

applications/solvers/multiphase/interIsoFoam/correctPhi.H

+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"

applications/solvers/multiphase/interIsoFoam/createFields.H

+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "createPhi.H"
+
+
+Info<< "Reading transportProperties\n" << endl;
+immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
+
+volScalarField& alpha1(mixture.alpha1());
+volScalarField& alpha2(mixture.alpha2());
+
+const dimensionedScalar& rho1 = mixture.rho1();
+const dimensionedScalar& rho2 = mixture.rho2();
+
+
+// Need to store rho for ddt(rho, U)
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT
+    ),
+    alpha1*rho1 + alpha2*rho2
+);
+rho.oldTime();
+
+
+// Mass flux
+surfaceScalarField rhoPhi
+(
+    IOobject
+    (
+        "rhoPhi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    fvc::interpolate(rho)*phi
+);
+
+
+// Construct incompressible turbulence model
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, mixture)
+);
+
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    p_rgh + rho*gh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell
+(
+    p,
+    p_rgh,
+    pimple.dict(),
+    pRefCell,
+    pRefValue
+);
+
+if (p_rgh.needReference())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+    p_rgh = p - rho*gh;
+}
+
+mesh.setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(alpha1.name());
+
+#include "createMRF.H"
+#include "createIsoAdvection.H"

applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H

+// Defining isoAdvection
+isoAdvection advector(alpha1, phi, U);
+
+bool frozenFlow = pimple.dict().lookupOrDefault<bool>("frozenFlow", false);
+if (frozenFlow)
+{
+    Info<< "Employing frozen-flow assumption: "
+        << "pressure-velocity system will not be updated"
+        << endl;
+}

applications/solvers/multiphase/interIsoFoam/interIsoFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016 DHI
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    interIsoFoam
+
+Group
+    grpMultiphaseSolvers
+
+Description
+    Solver derived from interFoam for 2 incompressible, isothermal immiscible
+    fluids using the iso-advector phase-fraction based interface capturing
+    approach.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+    For a two-fluid approach see twoPhaseEulerFoam.
+
+    Reference:
+        \verbatim
+            Roenby, J., Bredmose, H. and Jasak, H. (2016).
+            A computational method for sharp interface advection
+            Royal Society Open Science, 3
+            doi 10.1098/rsos.160405
+        \endverbatim
+
+    isoAdvector code supplied by Johan Roenby, DHI (2016)
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoAdvection.H"
+#include "fvCFD.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "correctPhi.H"
+
+    turbulence->validate();
+
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
+            mixture.correct();
+
+            if (frozenFlow)
+            {
+                continue;
+            }
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interIsoFoam/pEqn.H

+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+    MRF.makeRelative(phiHbyA);
+    adjustPhi(phiHbyA, U, p_rgh);
+
+    surfaceScalarField phig
+    (
+        (
+            mixture.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+//          - ghf*(fvc::grad(rho) & mesh.Sf())*rAUf
+    );
+
+    phiHbyA += phig;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+}

applications/solvers/multiphase/interIsoFoam/setDeltaT.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    setDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum courant Number.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    scalar maxDeltaTFact =
+        min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
+
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            deltaTFact*runTime.deltaTValue(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //

applications/utilities/preProcessing/setAlphaField/Make/files

+setAlphaField.C
+
+EXE = $(FOAM_APPBIN)/setAlphaField

applications/utilities/preProcessing/setAlphaField/Make/options

+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling

applications/utilities/preProcessing/setAlphaField/setAlphaField.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016-2017 DHI
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    setAlphaField
+
+Description
+    Uses isoCutCell to create a volume fraction field from either a cylinder,
+    a sphere or a plane.
+
+    Original code supplied by Johan Roenby, DHI (2016)
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "isoCutFace.H"
+#include "isoCutCell.H"
+#include "Enum.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+class shapeSelector
+{
+    public:
+
+        enum class shapeType
+        {
+            PLANE,
+            SPHERE,
+            CYLINDER,
+            SIN
+        };
+
+    static const Foam::Enum<shapeType> shapeTypeNames;
+};
+
+const Foam::Enum<shapeSelector::shapeType> shapeSelector::shapeTypeNames
+{
+    { shapeSelector::shapeType::PLANE, "plane" },
+    { shapeSelector::shapeType::SPHERE, "sphere" },
+    { shapeSelector::shapeType::CYLINDER, "cylinder" },
+    { shapeSelector::shapeType::SIN, "sin" },
+};
+
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    Info<< "Reading setAlphaFieldDict\n" << endl;
+
+    IOdictionary dict
+    (
+        IOobject
+        (
+            "setAlphaFieldDict",
+            runTime.system(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    const shapeSelector::shapeType surfType
+    (
+        shapeSelector::shapeTypeNames.read(dict.lookup("type"))
+    );
+    const vector centre(dict.lookup("centre"));
+    const word fieldName(dict.lookup("field"));
+
+    Info<< "Reading field " << fieldName << "\n" << endl;
+    volScalarField alpha1
+    (
+        IOobject
+        (
+            fieldName,
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    scalar f0 = 0.0;
+    scalarField f(mesh.points().size());
+
+    Info<< "Processing type '" << shapeSelector::shapeTypeNames[surfType]
+        << "'" << endl;
+
+    switch (surfType)
+    {
+        case shapeSelector::shapeType::PLANE:
+        {
+            const vector direction(dict.lookup("direction"));
+
+            f = -(mesh.points() - centre) & (direction/mag(direction));
+            f0 = 0.0;
+            break;
+        }
+        case shapeSelector::shapeType::SPHERE:
+        {
+            const scalar radius(readScalar(dict.lookup("radius")));
+
+            f = -mag(mesh.points() - centre);
+            f0 = -radius;
+            break;
+        }
+        case shapeSelector::shapeType::CYLINDER:
+        {
+            const scalar radius(readScalar(dict.lookup("radius")));
+            const vector direction(dict.lookup("direction"));
+
+            f = -sqrt
+            (
+                sqr(mag(mesh.points() - centre))
+              - sqr(mag((mesh.points() - centre) & direction))
+            );
+            f0 = -radius;
+            break;
+        }
+        case shapeSelector::shapeType::SIN:
+        {
+            const scalar period(readScalar(dict.lookup("period")));
+            const scalar amplitude(readScalar(dict.lookup("amplitude")));
+            const vector up(dict.lookup("up"));
+            const vector direction(dict.lookup("direction"));
+
+            const scalarField xx
+            (
+                (mesh.points() - centre) & direction/mag(direction)
+            );
+            const scalarField zz((mesh.points() - centre) & up/mag(up));
+
+            f = amplitude*Foam::sin(2*mathematical::pi*xx/period) - zz;
+            f0 = 0;
+            break;
+        }
+    }
+
+
+    // Define function on mesh points and isovalue
+
+    // Calculating alpha1 volScalarField from f = f0 isosurface
+    isoCutCell icc(mesh, f);
+    icc.volumeOfFluid(alpha1, f0);
+
+    // Writing volScalarField alpha1
+    ISstream::defaultPrecision(18);
+    alpha1.write();
+
+    Info<< nl << "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+        << "  Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
+        << nl << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+field           alpha.water;
+type            cylinder;
+radius          0.25;
+direction       (0 1 0);
+centre          (0.5 0 0.5);
+
+// ************************************************************************* //

src/finiteVolume/Make/files

@@ -238,6 +238,9 @@ fvMatrices/fvMatrices.C
 fvMatrices/fvScalarMatrix/fvScalarMatrix.C
 fvMatrices/solvers/MULES/MULES.C
 fvMatrices/solvers/MULES/CMULES.C
+fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C
+fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
+fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
 fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
 
 interpolation = interpolation/interpolation

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016-2017 DHI
+-------------------------------------------------------------------------------
+
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::isoAdvection
+
+Description
+    Calculates the new VOF (alpha) field after time step dt given the initial
+    VOF field and a velocity field U and face fluxes phi. The fluid transport
+    calculation is based on an idea of using isosurfaces to estimate the
+    internal distribution of fluid in cells and advecting such isosurfaces
+    across the mesh faces with the velocity field interpolated to the
+    isosurfaces.
+
+    Reference:
+        \verbatim
+            Roenby, J., Bredmose, H. and Jasak, H. (2016).
+            A computational method for sharp interface advection
+            Royal Society Open Science, 3
+            doi 10.1098/rsos.160405
+        \endverbatim
+
+    Original code supplied by Johan Roenby, DHI (2016)
+
+SourceFiles
+    isoAdvection.C
+    isoAdvectionTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoAdvection_H
+#define isoAdvection_H
+
+#include "fvMesh.H"
+#include "volFieldsFwd.H"
+#include "surfaceFields.H"
+#include "className.H"
+#include "isoCutCell.H"
+#include "isoCutFace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class isoAdvection Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoAdvection
+{
+    // Private data types
+
+        typedef DynamicList<label> DynamicLabelList;
+        typedef DynamicList<scalar> DynamicScalarList;
+        typedef DynamicList<vector> DynamicVectorList;
+        typedef DynamicList<point> DynamicPointList;
+
+
+    // Private data
+
+        //- Reference to mesh
+        const fvMesh& mesh_;
+
+        //- Dictionary for isoAdvection controls
+        const dictionary dict_;
+
+        //- VOF field
+        volScalarField& alpha1_;
+
+        //- Often used reference to alpha1 internal field
+        scalarField& alpha1In_;
+
+        //- Reference to flux field
+        const surfaceScalarField& phi_;
+
+        //- Reference to velocity field
+        const volVectorField& U_;
+
+        //- Face volumetric water transport
+        surfaceScalarField dVf_;
+
+        //- Time spent performing interface advection
+        scalar advectionTime_;
+
+
+        // Point interpolation data
+
+            //- VOF field interpolated to mesh points
+            scalarField ap_;
+
+
+        // Switches and tolerances. Tolerances need to go into toleranceSwitches
+
+            //- Number of alpha bounding steps
+            label nAlphaBounds_;
+
+            //- Tolerance for search of isoFace giving specified VOF value
+            scalar vof2IsoTol_;
+
+            //- Tolerance for marking of surface cells:
+            //  Those with surfCellTol_ < alpha1 < 1 - surfCellTol_
+            scalar surfCellTol_;
+
+            //- Switch controlling whether to use isoface normals for interface
+            //  orientation (default corresponding to false) to base it on
+            //  a smoothed gradient of alpha calculation (giving better results
+            //  on tri on tet meshes).
+            bool gradAlphaBasedNormal_;
+
+            //- Print isofaces in a <case>/isoFaces/isoFaces_#N.vtk files.
+            //  Intended for debugging
+            bool writeIsoFacesToFile_;
+
+        // Cell and face cutting
+
+            //- List of surface cells
+            DynamicLabelList surfCells_;
+
+            //- Cell cutting object
+            isoCutCell isoCutCell_;
+
+            //- Face cutting object
+            isoCutFace isoCutFace_;
+
+            //- Bool list for cells that have been touched by the bounding step
+            boolList cellIsBounded_;
+
+            //- True for all surface cells and their neighbours
+            boolList checkBounding_;
+
+            //- Storage for boundary faces downwind to a surface cell
+            DynamicLabelList bsFaces_;
+
+            //- Storage for boundary surface iso face centre
+            DynamicVectorList bsx0_;
+
+            //- Storage for boundary surface iso face normal
+            DynamicVectorList bsn0_;
+
+            //- Storage for boundary surface iso face speed
+            DynamicScalarList bsUn0_;
+
+            //- Storage for boundary surface iso value
+            DynamicScalarList bsf0_;
+
+
+        // Additional data for parallel runs
+
+            //- List of processor patch labels
+            DynamicLabelList procPatchLabels_;
+
+            //- For each patch if it is a processor patch this is a list of the
+            //  face labels on this patch that are downwind to a surface cell.
+            //  For non-processor patches the list will be empty.
+            List<DynamicLabelList> surfaceCellFacesOnProcPatches_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        isoAdvection(const isoAdvection&);
+
+        //- Disallow default bitwise copy assignment
+        void operator=(const isoAdvection&);
+
+
+        // Advection functions
+
+            //- For each face calculate volumetric face transport during dt
+            void timeIntegratedFlux();
+
+            //- Calculate volumetric face transport during dt given the isoFace
+            //  data provided as input for face facei
+            scalar timeIntegratedFaceFlux
+            (
+                const label facei,
+                const vector& x0,
+                const vector& n0,
+                const scalar Un0,
+                const scalar f0,
+                const scalar dt,
+                const scalar phi,
+                const scalar magSf
+            );
+
+            //- For a given cell return labels of faces fluxing out of this cell
+            //  (based on sign of phi)
+            void setDownwindFaces
+            (
+                const label celli,
+                DynamicLabelList& downwindFaces
+            ) const;
+
+            // Limit fluxes
+            void limitFluxes();
+
+            // Bound fluxes
+            void boundFromAbove
+            (
+                const scalarField& alpha1,
+                surfaceScalarField& dVfcorrected,
+                DynamicLabelList& correctedFaces
+            );
+
+            //- Given the face volume transport dVf calculates the total volume
+            //  leaving a given cell. Note: cannot use dVf member because
+            //  netFlux is called also for corrected dVf
+            scalar netFlux
+            (
+                const surfaceScalarField& dVf,
+                const label celli
+            ) const;
+
+            //- Determine if a cell is a surface cell
+            bool isASurfaceCell(const label celli) const
+            {
+                return
+                (
+                    surfCellTol_ < alpha1In_[celli]
+                 && alpha1In_[celli] < 1 - surfCellTol_
+                );
+            }
+
+            //- Clear out isoFace data
+            void clearIsoFaceData()
+            {
+                surfCells_.clear();
+                bsFaces_.clear();
+                bsx0_.clear();
+                bsn0_.clear();
+                bsUn0_.clear();
+                bsf0_.clear();
+            }
+
+        // Face value functions needed for random face access where the face
+        // can be either internal or boundary face
+
+            //- Return face value for a given Geometric surface field
+            template<typename Type>
+            Type faceValue
+            (
+                const GeometricField<Type, fvsPatchField, surfaceMesh>& f,
+                const label facei
+            ) const;
+
+            //- Set face value for a given Geometric surface field
+            template<typename Type>
+            void setFaceValue
+            (
+                GeometricField<Type, fvsPatchField, surfaceMesh>& f,
+                const label facei,
+                const Type& value
+            ) const;
+
+
+        // Parallel run handling functions
+
+            //- Synchronize dVf across processor boundaries using upwind value
+            void syncProcPatches
+            (
+                surfaceScalarField& dVf,
+                const surfaceScalarField& phi
+            );
+
+            //- Check if the face is on processor patch and append it to the
+            //  list of surface cell faces on processor patches
+            void checkIfOnProcPatch(const label facei);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("isoAdvection");
+
+    //- Constructors
+
+        //- Construct given alpha, phi and velocity field. Note: phi should be
+        // divergence free up to a sufficient tolerance
+        isoAdvection
+        (
+            volScalarField& alpha1,
+            const surfaceScalarField& phi,
+            const volVectorField& U
+        );
+
+
+    //- Destructor
+    virtual ~isoAdvection()
+    {}
+
+
+    // Member functions
+
+        //- Advect the free surface. Updates alpha field, taking into account
+        //  multiple calls within a single time step.
+        void advect();
+
+        //- Apply the bounding based on user inputs
+        void applyBruteForceBounding();
+
+        // Access functions
+
+            //- Return alpha field
+            const volScalarField& alpha() const
+            {
+                return alpha1_;
+            }
+
+            //- Return the controls dictionary
+            const dictionary& dict() const
+            {
+                return dict_;
+            }
+
+            //- Return cellSet of surface cells
+            void writeSurfaceCells() const;
+
+            //- Return cellSet of bounded cells
+            void writeBoundedCells() const;
+
+            //- Return mass flux
+            tmp<surfaceScalarField> getRhoPhi
+            (
+                const dimensionedScalar rho1,
+                const dimensionedScalar rho2
+            ) const
+            {
+                return tmp<surfaceScalarField>
+                (
+                    new surfaceScalarField
+                    (
+                        "rhoPhi",
+                        (rho1 - rho2)*dVf_/mesh_.time().deltaT() + rho2*phi_
+                    )
+                );
+            }
+
+            scalar advectionTime() const
+            {
+                return advectionTime_;
+            }
+
+            //- Write isoface points to .obj file
+            void writeIsoFaces
+            (
+                const DynamicList<List<point>>& isoFacePts
+            ) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "isoAdvectionTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //