diff --git a/applications/solvers/multiphase/interIsoFoam/Make/options b/applications/solvers/multiphase/interIsoFoam/Make/options
index 21e9667b83..6d6417f7a0 100644
--- a/applications/solvers/multiphase/interIsoFoam/Make/options
+++ b/applications/solvers/multiphase/interIsoFoam/Make/options
@@ -7,6 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
@@ -15,6 +16,7 @@ EXE_LIBS = \
     -lturbulenceModels \
     -lincompressibleTurbulenceModels \
     -lfiniteVolume \
+    -ldynamicFvMesh \
     -lfvOptions \
     -lmeshTools \
     -lsampling \
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaControls.H b/applications/solvers/multiphase/interIsoFoam/alphaControls.H
index db77d94af4..ebe73c3c2c 100644
--- a/applications/solvers/multiphase/interIsoFoam/alphaControls.H
+++ b/applications/solvers/multiphase/interIsoFoam/alphaControls.H
@@ -1,3 +1,29 @@
 const dictionary& alphaControls = mesh.solverDict(alpha1.name());
 
+//label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
+
 label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+
+/*
+bool MULESCorr(alphaControls.lookupOrDefault("MULESCorr", false));
+
+// Apply the compression correction from the previous iteration
+// Improves efficiency for steady-simulations but can only be applied
+// once the alpha field is reasonably steady, i.e. fully developed
+bool alphaApplyPrevCorr
+(
+    alphaControls.lookupOrDefault("alphaApplyPrevCorr", false)
+);
+
+// Isotropic compression coefficient
+scalar icAlpha
+(
+    alphaControls.lookupOrDefault<scalar>("icAlpha", 0)
+);
+
+// Shear compression coefficient
+scalar scAlpha
+(
+    alphaControls.lookupOrDefault<scalar>("scAlpha", 0)
+);
+*/
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H b/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H
index a084155c8a..24f6e48a22 100644
--- a/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H
+++ b/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H
index cefabfda73..1bda9eef18 100644
--- a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H
@@ -1,35 +1,50 @@
-// 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();
+    {
+        // Note for AMR implementation:
+        // At this point we have just entered the new time step,
+        // the mesh has been refined and the alpha, phi and U contain
+        // the field values at the beginning of this time step mapped
+        // to the new mesh.
 
-    // 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());
-}
+        // The above assumes that we are in firstIter() of the outer
+        // corrector loop. If we are in any subsequent iter of this loop
+        // the alpha1, U and phi will be overwritten with the new time step
+        // values but still on the same mesh.
 
-// Update alpha1
-advector.advect();
 
-// Update rhoPhi
-rhoPhi = advector.getRhoPhi(rho1, rho2);
+        if (pimple.firstIter())
+        {
+            // Note: This assumes moveMeshOuterCorrectors is false
+            alpha10 = alpha1;
+            U0 = U;
+            phi0 = phi;
+            advector.advect();
 
-alpha2 = 1.0 - alpha1;
+            #include "rhofs.H"
+            rhoPhi = advector.getRhoPhi(rho1f, rho2f);
+        }
+        else
+        {
+            alpha1 = alpha10;
+            // Temporarily setting U and phi to average of old and new value
+            // Note: Illegal additions if mesh is topoChanging
+            // (e.g. if moveMeshOuterCorrectors and AMR)
+            U = 0.5*U0 + 0.5*U;
+            phi = 0.5*phi0 + 0.5*phi;
+            isoAdvection advector(alpha1, phi, U);
+            advector.advect();
+            #include "rhofs.H"
+            rhoPhi = advector.getRhoPhi(rho1f, rho2f);
+            // Resetting U and phi to the new value
+            U = 2.0*U - U0;
+            phi = 2.0*phi - phi0;
+        }
 
-if (!pimple.firstIter())
-{
-    // Restoring new U and phi values temporarily saved in prevIter() above
-    U = U.prevIter();
-    phi = phi.prevIter();
-}
+        alpha2 = 1.0 - alpha1;
+        mixture.correct();
 
-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;
+    }
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+        << endl;
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H
index 8fd965e312..c53d6e106b 100644
--- a/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H
@@ -15,6 +15,12 @@ if (nAlphaSubCycles > 1)
 
     tmp<volScalarField> trSubDeltaT;
 
+    if (LTS)
+    {
+        trSubDeltaT =
+            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+    }
+
     for
     (
         subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaSuSp.H b/applications/solvers/multiphase/interIsoFoam/alphaSuSp.H
new file mode 100644
index 0000000000..0b8e05e187
--- /dev/null
+++ b/applications/solvers/multiphase/interIsoFoam/alphaSuSp.H
@@ -0,0 +1,3 @@
+zeroField Su;
+zeroField Sp;
+zeroField divU;
diff --git a/applications/solvers/multiphase/interIsoFoam/correctPhi.H b/applications/solvers/multiphase/interIsoFoam/correctPhi.H
index 9afcd58a66..fcb5020587 100644
--- a/applications/solvers/multiphase/interIsoFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interIsoFoam/correctPhi.H
@@ -3,7 +3,7 @@ CorrectPhi
     U,
     phi,
     p_rgh,
-    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
     geometricZeroField(),
     pimple
 );
diff --git a/applications/solvers/multiphase/interIsoFoam/createFields.H b/applications/solvers/multiphase/interIsoFoam/createFields.H
index 8c3a05e3fa..02269e379a 100644
--- a/applications/solvers/multiphase/interIsoFoam/createFields.H
+++ b/applications/solvers/multiphase/interIsoFoam/createFields.H
@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
 (
@@ -119,5 +121,67 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
+/*
+// MULES compressed flux is registered in case scalarTransport FO needs it.
+surfaceScalarField alphaPhiUn
+(
+    IOobject
+    (
+        "alphaPhiUn",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(phi.dimensions(), Zero)
+);
+*/
+
 #include "createMRF.H"
-#include "createIsoAdvection.H"
+#include "createFvOptions.H"
+
+//Auxiliary fields when nOuterCorrectors > 1
+surfaceScalarField phi0
+(
+    IOobject
+    (
+        "phi0",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    phi
+);
+
+
+volVectorField U0
+(
+    IOobject
+    (
+        "U0",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    U
+);
+
+
+volScalarField alpha10
+(
+    IOobject
+    (
+        "alpha10",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    alpha1
+);
+
+
+isoAdvection advector(alpha1, phi, U);
diff --git a/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H b/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H
deleted file mode 100644
index 3360cdae65..0000000000
--- a/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H
+++ /dev/null
@@ -1,2 +0,0 @@
-// Defining isoAdvection
-isoAdvection advector(alpha1, phi, U);
diff --git a/applications/solvers/multiphase/interIsoFoam/initCorrectPhi.H b/applications/solvers/multiphase/interIsoFoam/initCorrectPhi.H
new file mode 100644
index 0000000000..03eb7fd7fd
--- /dev/null
+++ b/applications/solvers/multiphase/interIsoFoam/initCorrectPhi.H
@@ -0,0 +1,34 @@
+tmp<volScalarField> rAU;
+
+if (correctPhi)
+{
+    rAU = new volScalarField
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rAU", dimTime/dimDensity, 1)
+    );
+
+    #include "correctPhi.H"
+}
+else
+{
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+        geometricZeroField(),
+        pimple
+    );
+
+    #include "continuityErrs.H"
+}
diff --git a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C
index 07a0cf1e6d..024e502633 100644
--- a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C
+++ b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
                 isoAdvector | Copyright (C) 2016 DHI
+              Modified work | Copyright (C) 2018 Johan Roenby
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,36 +32,35 @@ Group
 
 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.
+    fluids using the isoAdvector phase-fraction based interface capturing
+    approach, with optional mesh motion and mesh topology changes including
+    adaptive re-meshing.
 
     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
+    \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)
+    isoAdvector code supplied by Johan Roenby, STROMNING (2018)
 
 \*---------------------------------------------------------------------------*/
 
-#include "isoAdvection.H"
 #include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "isoAdvection.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
 #include "immiscibleIncompressibleTwoPhaseMixture.H"
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
+#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,31 +71,39 @@ int main(int argc, char *argv[])
     #include "addCheckCaseOptions.H"
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createTimeControls.H"
+    #include "createDynamicFvMesh.H"
     #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
     #include "createFields.H"
-    #include "createFvOptions.H"
-    #include "correctPhi.H"
+//    #include "createAlphaFluxes.H"
+    #include "initCorrectPhi.H"
+    #include "createUfIfPresent.H"
 
     turbulence->validate();
 
-    #include "readTimeControls.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
+    if (!LTS)
+    {
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "readDyMControls.H"
 
-        #include "CourantNo.H"
-        #include "alphaCourantNo.H"
-        #include "setDeltaT.H"
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "CourantNo.H"
+            #include "alphaCourantNo.H"
+            #include "setDeltaT.H"
+        }
 
         runTime++;
 
@@ -104,6 +112,45 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    // Do not apply previous time-step mesh compression flux
+                    // if the mesh topology changed
+                    if (mesh.topoChanging())
+                    {
+//                        talphaPhi1Corr0.clear();
+                    }
+
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+
+                    MRF.update();
+
+                    if (correctPhi)
+                    {
+                        // Calculate absolute flux
+                        // from the mapped surface velocity
+                        phi = mesh.Sf() & Uf();
+
+                        #include "correctPhi.H"
+
+                        // Make the flux relative to the mesh motion
+                        fvc::makeRelative(phi, U);
+
+                        mixture.correct();
+                    }
+
+                    if (checkMeshCourantNo)
+                    {
+                        #include "meshCourantNo.H"
+                    }
+                }
+            }
+
             #include "alphaControls.H"
             #include "alphaEqnSubCycle.H"
 
diff --git a/applications/solvers/multiphase/interIsoFoam/pEqn.H b/applications/solvers/multiphase/interIsoFoam/pEqn.H
index 17c3ae9748..e443f26984 100644
--- a/applications/solvers/multiphase/interIsoFoam/pEqn.H
+++ b/applications/solvers/multiphase/interIsoFoam/pEqn.H
@@ -1,17 +1,29 @@
 {
-    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+    if (correctPhi)
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        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.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
     );
     MRF.makeRelative(phiHbyA);
-    adjustPhi(phiHbyA, U, p_rgh);
+
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
 
     surfaceScalarField phig
     (
@@ -19,7 +31,6 @@
             mixture.surfaceTensionForce()
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf()
-//          - ghf*(fvc::grad(rho) & mesh.Sf())*rAUf
     );
 
     phiHbyA += phig;
@@ -44,7 +55,7 @@
 
             p_rgh.relax();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
@@ -52,6 +63,12 @@
 
     #include "continuityErrs.H"
 
+    // Correct Uf if the mesh is moving
+    fvc::correctUf(Uf, U, phi);
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
     p == p_rgh + rho*gh;
 
     if (p_rgh.needReference())
@@ -64,4 +81,9 @@
         );
         p_rgh = p - rho*gh;
     }
+
+    if (!correctPhi)
+    {
+        rAU.clear();
+    }
 }
diff --git a/applications/solvers/multiphase/interIsoFoam/rhofs.H b/applications/solvers/multiphase/interIsoFoam/rhofs.H
new file mode 100644
index 0000000000..5949bf1e0a
--- /dev/null
+++ b/applications/solvers/multiphase/interIsoFoam/rhofs.H
@@ -0,0 +1,2 @@
+const dimensionedScalar& rho1f(rho1);
+const dimensionedScalar& rho2f(rho2);
diff --git a/applications/solvers/multiphase/interIsoFoam/setDeltaT.H b/applications/solvers/multiphase/interIsoFoam/setDeltaT.H
index 9cc860c032..924d24c8ea 100644
--- a/applications/solvers/multiphase/interIsoFoam/setDeltaT.H
+++ b/applications/solvers/multiphase/interIsoFoam/setDeltaT.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/solvers/multiphase/interIsoFoam/setRDeltaT.H b/applications/solvers/multiphase/interIsoFoam/setRDeltaT.H
new file mode 100644
index 0000000000..26b237e795
--- /dev/null
+++ b/applications/solvers/multiphase/interIsoFoam/setRDeltaT.H
@@ -0,0 +1,136 @@
+{
+    volScalarField& rDeltaT = trDeltaT.ref();
+
+    const dictionary& pimpleDict = pimple.dict();
+
+    scalar maxCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
+    );
+
+    scalar maxAlphaCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
+    );
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    label nAlphaSpreadIter
+    (
+        pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
+    );
+
+    scalar alphaSpreadDiff
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
+    );
+
+    scalar alphaSpreadMax
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
+    );
+
+    scalar alphaSpreadMin
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
+    );
+
+    label nAlphaSweepIter
+    (
+        pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
+    );
+
+    scalar rDeltaTDampingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+    );
+
+    scalar maxDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    );
+
+    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.ref() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(mag(rhoPhi))()()
+       /((2*maxCo)*mesh.V()*rho())
+    );
+
+    if (maxAlphaCo < maxCo)
+    {
+        // Further limit the reciprocal time-step
+        // in the vicinity of the interface
+
+        volScalarField alpha1Bar(fvc::average(alpha1));
+
+        rDeltaT.ref() = max
+        (
+            rDeltaT(),
+            pos0(alpha1Bar() - alphaSpreadMin)
+           *pos0(alphaSpreadMax - alpha1Bar())
+           *fvc::surfaceSum(mag(phi))()()
+           /((2*maxAlphaCo)*mesh.V())
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    if (nAlphaSpreadIter > 0)
+    {
+        fvc::spread
+        (
+            rDeltaT,
+            alpha1,
+            nAlphaSpreadIter,
+            alphaSpreadDiff,
+            alphaSpreadMax,
+            alphaSpreadMin
+        );
+    }
+
+    if (nAlphaSweepIter > 0)
+    {
+        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
+    }
+
+    Info<< "Smoothed flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+
+        Info<< "Damped flow time scale min/max = "
+            << gMin(1/rDeltaT.primitiveField())
+            << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    }
+}
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
index 3382092993..1d4062a998 100644
--- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
@@ -156,7 +156,6 @@ void Foam::isoAdvection::timeIntegratedFlux()
     // Clear out the data for re-use and reset list containing information
     // whether cells could possibly need bounding
     clearIsoFaceData();
-    checkBounding_ = false;
 
     // Get necessary references
     const scalarField& phiIn = phi_.primitiveField();
@@ -164,201 +163,172 @@ void Foam::isoAdvection::timeIntegratedFlux()
     scalarField& dVfIn = dVf_.primitiveFieldRef();
 
     // Get necessary mesh data
-    const labelListList& cellPoints = mesh_.cellPoints();
-    const labelListList& cellCells = mesh_.cellCells();
     const cellList& cellFaces = mesh_.cells();
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
-    const vectorField& cellCentres = mesh_.cellCentres();
-    const pointField& points = mesh_.points();
+    const labelListList& cellCells = mesh_.cellCells();
 
-    // Storage for isoFace points. Only used if writeIsoFacesToFile_
-    DynamicList<List<point>> isoFacePts;
-
-    // Interpolating alpha1 cell centre values to mesh points (vertices)
-    ap_ = volPointInterpolation::New(mesh_).interpolate(alpha1_);
-
-    vectorField gradAlpha(mesh_.nPoints(), vector::zero);
+    // Calculate alpha vertex values, ap_, or cell normals (used to get
+    // interface-vertex distance function if gradAlphaBasedNormal_)
+    volVectorField cellNormals("cellN", fvc::grad(alpha1_));
     if (gradAlphaBasedNormal_)
     {
-        // Calculate gradient of alpha1 and interpolate to vertices
-        volVectorField gradA("gradA", fvc::grad(alpha1_));
-        gradAlpha = volPointInterpolation::New(mesh_).interpolate(gradA);
+        // Calculate gradient of alpha1 and normalise and smoothen it.
+        normaliseAndSmooth(cellNormals);
+    }
+    else
+    {
+        // Interpolating alpha1 cell centre values to mesh points (vertices)
+        ap_ = volPointInterpolation::New(mesh_).interpolate(alpha1_);
     }
 
-    // Loop through cells
+    // Storage for isoFace points. Only used if writeIsoFacesToFile_
+    DynamicList<List<point> > isoFacePts;
+
+    // Loop through all cells
     forAll(alpha1In_, celli)
     {
-        if (isASurfaceCell(celli))
+        // If not a surface cell continue to next cell
+        if (!isASurfaceCell(celli)) continue;
+
+        // This is a surface cell, increment counter, append and mark cell
+        // Note: We also have the cellStatus below where the cell might not have
+        // an isoface. So maybe the counter and append should be put there.
+        nSurfaceCells++;
+        surfCells_.append(celli);
+        checkBounding_[celli] = true;
+
+        DebugInfo
+            << "\n------------ Cell " << celli << " with alpha1 = "
+            << alpha1In_[celli] << " and 1-alpha1 = "
+            << 1.0 - alpha1In_[celli] << " ------------"
+            << endl;
+
+        if (gradAlphaBasedNormal_)
         {
-            // This is a surface cell, increment counter, append and mark cell
-            nSurfaceCells++;
-            surfCells_.append(celli);
-            checkBounding_.set(celli);
+            vectorField& cellNormalsIn = cellNormals.primitiveFieldRef();
+            setCellVertexValues(celli, cellNormalsIn);
+        }
 
-            DebugInfo
-                << "\n------------ Cell " << celli << " with alpha1 = "
-                << alpha1In_[celli] << " and 1-alpha1 = "
-                << 1.0 - alpha1In_[celli] << " ------------"
-                << endl;
+        // Calculate isoFace centre x0, normal n0 at time t
 
-            // Calculate isoFace centre x0, normal n0 at time t
-            label maxIter = 100; // NOTE: make it a debug switch
+        // Calculate cell status (-1: cell is fully below the isosurface, 0:
+        // cell is cut, 1: cell is fully above the isosurface)
+        label maxIter = 100; // NOTE: make it a debug switch
+        label cellStatus = isoCutCell_.vofCutCell
+        (
+            celli,
+            alpha1In_[celli],
+            isoFaceTol_,
+            maxIter
+        );
 
-            const labelList& cp = cellPoints[celli];
-            scalarField ap_org(cp.size(), 0);
-            if (gradAlphaBasedNormal_)
+        // If cell is not cut move on to next cell
+        if (cellStatus != 0) continue;
+
+        // If cell is cut calculate isoface unit normal
+        const scalar f0(isoCutCell_.isoValue());
+        const point& x0(isoCutCell_.isoFaceCentre());
+        vector n0(isoCutCell_.isoFaceArea());
+        n0 /= (mag(n0));
+
+        if (writeIsoFacesToFile_ && mesh_.time().writeTime())
+        {
+            isoFacePts.append(isoCutCell_.isoFacePoints());
+        }
+
+        // Get the speed of the isoface by interpolating velocity and
+        // dotting it with isoface unit normal
+        const scalar Un0 = UInterp.interpolate(x0, celli) & n0;
+
+        DebugInfo
+            << "calcIsoFace gives initial surface: \nx0 = " << x0
+            << ", \nn0 = " << n0 << ", \nf0 = " << f0 << ", \nUn0 = "
+            << Un0 << endl;
+
+        // Estimate time integrated flux through each downwind face
+        // Note: looping over all cell faces - in reduced-D, some of
+        //       these faces will be on empty patches
+        const cell& celliFaces = cellFaces[celli];
+        forAll(celliFaces, fi)
+        {
+            const label facei = celliFaces[fi];
+
+            if (mesh_.isInternalFace(facei))
             {
-                // Calculating smoothed alpha gradient in surface cell in order
-                // to use it as the isoface orientation.
-                vector smoothedGradA = vector::zero;
-                const point& cellCentre = cellCentres[celli];
-                scalar wSum = 0;
-                forAll(cp, pointI)
-                {
-                    point vertex = points[cp[pointI]];
-                    scalar w = 1.0/mag(vertex - cellCentre);
-                    wSum += w;
-                    smoothedGradA += w*gradAlpha[cp[pointI]];
-                }
-                smoothedGradA /= wSum;
+                bool isDownwindFace = false;
+                label otherCell = -1;
 
-                // Temporarily overwrite the interpolated vertex alpha values in
-                // ap_ with the vertex-cell centre distance along smoothedGradA.
-                forAll(ap_org, vi)
+                if (celli == own[facei])
                 {
-                    ap_org[vi] = ap_[cp[vi]];
-                    const point& vertex = points[cp[vi]];
-                    ap_[cp[vi]] =
+                    if (phiIn[facei] > 10*SMALL)
+                    {
+                        isDownwindFace = true;
+                    }
+
+                    otherCell = nei[facei];
+                }
+                else
+                {
+                    if (phiIn[facei] < -10*SMALL)
+                    {
+                        isDownwindFace = true;
+                    }
+
+                    otherCell = own[facei];
+                }
+
+                if (isDownwindFace)
+                {
+                    dVfIn[facei] = isoCutFace_.timeIntegratedFaceFlux
                     (
-                        (vertex - cellCentre)
-                      & (smoothedGradA/mag(smoothedGradA))
+                        facei,
+                        x0,
+                        n0,
+                        Un0,
+                        f0,
+                        dt,
+                        phiIn[facei],
+                        magSfIn[facei]
                     );
                 }
-            }
 
-            // Calculate cell status (-1: cell is fully below the isosurface, 0:
-            // cell is cut, 1: cell is fully above the isosurface)
-            label cellStatus = isoCutCell_.vofCutCell
-            (
-                celli,
-                alpha1In_[celli],
-                isoFaceTol_,
-                maxIter
-            );
+                // We want to check bounding of neighbour cells to
+                // surface cells as well:
+                checkBounding_[otherCell] = true;
 
-            if (gradAlphaBasedNormal_)
-            {
-                // Restoring ap_ by putting the original values  back into it.
-                forAll(ap_org, vi)
+                // Also check neighbours of neighbours.
+                // Note: consider making it a run time selectable
+                // extension level (easily done with recursion):
+                // 0 - only neighbours
+                // 1 - neighbours of neighbours
+                // 2 - ...
+                // Note: We will like all point neighbours to interface cells to
+                // be checked. Especially if the interface leaves a cell during
+                // a time step, it may enter a point neighbour which should also
+                // be treated like a surface cell. Its interface normal should 
+                // somehow be inherrited from its upwind cells from which it 
+                // receives the interface.
+                const labelList& nNeighbourCells = cellCells[otherCell];
+                forAll(nNeighbourCells, ni)
                 {
-                    ap_[cp[vi]] = ap_org[vi];
+                    checkBounding_[nNeighbourCells[ni]] = true;
                 }
             }
-
-            // Cell is cut
-            if (cellStatus == 0)
+            else
             {
-                const scalar f0 = isoCutCell_.isoValue();
-                const point& x0 = isoCutCell_.isoFaceCentre();
-                vector n0 = isoCutCell_.isoFaceArea();
-                n0 /= (mag(n0));
+                bsFaces_.append(facei);
+                bsx0_.append(x0);
+                bsn0_.append(n0);
+                bsUn0_.append(Un0);
+                bsf0_.append(f0);
 
-                if (writeIsoFacesToFile_ && mesh_.time().writeTime())
-                {
-                    isoFacePts.append(isoCutCell_.isoFacePoints());
-                }
-
-                // Get the speed of the isoface by interpolating velocity and
-                // dotting it with isoface normal
-                const scalar Un0 = UInterp.interpolate(x0, celli) & n0;
-
-                DebugInfo
-                    << "calcIsoFace gives initial surface: \nx0 = " << x0
-                    << ", \nn0 = " << n0 << ", \nf0 = " << f0 << ", \nUn0 = "
-                    << Un0 << endl;
-
-                // Estimate time integrated flux through each downwind face
-                // Note: looping over all cell faces - in reduced-D, some of
-                //       these faces will be on empty patches
-                const cell& celliFaces = cellFaces[celli];
-                forAll(celliFaces, fi)
-                {
-                    const label facei = celliFaces[fi];
-
-                    if (mesh_.isInternalFace(facei))
-                    {
-                        bool isDownwindFace = false;
-                        label otherCell = -1;
-
-                        if (celli == own[facei])
-                        {
-                            if (phiIn[facei] > 10*SMALL)
-                            {
-                                isDownwindFace = true;
-                            }
-
-                            otherCell = nei[facei];
-                        }
-                        else
-                        {
-                            if (phiIn[facei] < -10*SMALL)
-                            {
-                                isDownwindFace = true;
-                            }
-
-                            otherCell = own[facei];
-                        }
-
-                        if (isDownwindFace)
-                        {
-                            dVfIn[facei] = timeIntegratedFaceFlux
-                            (
-                                facei,
-                                x0,
-                                n0,
-                                Un0,
-                                f0,
-                                dt,
-                                phiIn[facei],
-                                magSfIn[facei]
-                            );
-                        }
-
-                        // We want to check bounding of neighbour cells to
-                        // surface cells as well:
-                        checkBounding_.set(otherCell);
-
-                        // Also check neighbours of neighbours.
-                        // Note: consider making it a run time selectable
-                        // extension level (easily done with recursion):
-                        // 0 - only neighbours
-                        // 1 - neighbours of neighbours
-                        // 2 - ...
-                        const labelList& nNeighbourCells = cellCells[otherCell];
-                        checkBounding_.set(nNeighbourCells);
-                    }
-                    else
-                    {
-                        bsFaces_.append(facei);
-                        bsx0_.append(x0);
-                        bsn0_.append(n0);
-                        bsUn0_.append(Un0);
-                        bsf0_.append(f0);
-
-                        // Note: we must not check if the face is on the
-                        // processor patch here.
-                    }
-                }
+                // Note: we must not check if the face is on the
+                // processor patch here.
             }
         }
     }
 
-    if (writeIsoFacesToFile_ && mesh_.time().writeTime())
-    {
-        writeIsoFaces(isoFacePts);
-    }
-
     // Get references to boundary fields
     const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
     const surfaceScalarField::Boundary& phib = phi_.boundaryField();
@@ -383,7 +353,7 @@ void Foam::isoAdvection::timeIntegratedFlux()
             {
                 const scalar magSf = magSfb[patchi][patchFacei];
 
-                dVfb[patchi][patchFacei] = timeIntegratedFaceFlux
+                dVfb[patchi][patchFacei] = isoCutFace_.timeIntegratedFaceFlux
                 (
                     facei,
                     bsx0_[i],
@@ -402,156 +372,62 @@ void Foam::isoAdvection::timeIntegratedFlux()
         }
     }
 
+    // Synchronize processor patches
+    syncProcPatches(dVf_, phi_);
+
+    writeIsoFaces(isoFacePts);
+
     Info<< "Number of isoAdvector surface cells = "
         << returnReduce(nSurfaceCells, sumOp<label>()) << endl;
 }
 
 
-Foam::scalar Foam::isoAdvection::timeIntegratedFaceFlux
+void Foam::isoAdvection::setCellVertexValues
 (
-    const label facei,
-    const vector& x0,
-    const vector& n0,
-    const scalar Un0,
-    const scalar f0,
-    const scalar dt,
-    const scalar phi,
-    const scalar magSf
+    const label celli,
+    const vectorField& cellNormalsIn
 )
 {
-    // Treating rare cases where isoface normal is not calculated properly
-    if (mag(n0) < 0.5)
+    const labelListList& cellPoints = mesh_.cellPoints();
+    const vectorField& cellCentres = mesh_.cellCentres();
+    const pointField& points = mesh_.points();
+    const labelList& cp = cellPoints[celli];
+    const point& cellCentre = cellCentres[celli];
+    forAll(cp, vi)
     {
-        scalar alphaf = 0;
-        scalar waterInUpwindCell = 0;
-
-        if (phi > 10*SMALL || !mesh_.isInternalFace(facei))
-        {
-            const label upwindCell = mesh_.faceOwner()[facei];
-            alphaf = alpha1In_[upwindCell];
-            waterInUpwindCell = alphaf*mesh_.cellVolumes()[upwindCell];
-        }
-        else
-        {
-            const label upwindCell = mesh_.faceNeighbour()[facei];
-            alphaf = alpha1In_[upwindCell];
-            waterInUpwindCell = alphaf*mesh_.cellVolumes()[upwindCell];
-        }
-
-        if (debug)
-        {
-            WarningInFunction
-                << "mag(n0) = " << mag(n0)
-                << " so timeIntegratedFlux calculates dVf from upwind"
-                << " cell alpha value: " << alphaf << endl;
-        }
-
-        return min(alphaf*phi*dt, waterInUpwindCell);
+        const point& vertex = points[cp[vi]];
+        ap_[cp[vi]] = (vertex - cellCentre) & cellNormalsIn[celli];
     }
+}
 
 
-    // Find sorted list of times where the isoFace will arrive at face points
-    // given initial position x0 and velocity Un0*n0
+void Foam::isoAdvection::normaliseAndSmooth
+(
+    volVectorField& cellN
+)
+{
+    const labelListList& cellPoints = mesh_.cellPoints();
+    const vectorField& cellCentres = mesh_.cellCentres();
+    const pointField& points = mesh_.points();
 
-    // Get points for this face
-    const face& f = mesh_.faces()[facei];
-    const pointField fPts(f.points(mesh_.points()));
-    const label nPoints = fPts.size();
-
-    scalarField pTimes(fPts.size());
-    if (mag(Un0) > 10*SMALL) // Note: tolerances
+    vectorField& cellNIn = cellN.primitiveFieldRef();
+    cellNIn /= (mag(cellNIn) + SMALL);
+    vectorField vertexN(mesh_.nPoints(), vector::zero);
+    vertexN = volPointInterpolation::New(mesh_).interpolate(cellN);
+    vertexN /= (mag(vertexN) + SMALL);
+    // Interpolate vertex normals back to cells
+    forAll(cellNIn, celli)
     {
-        // Here we estimate time of arrival to the face points from their normal
-        // distance to the initial surface and the surface normal velocity
-
-        pTimes = ((fPts - x0) & n0)/Un0;
-
-        scalar dVf = 0;
-
-        // Check if pTimes changes direction more than twice when looping face
-        label nShifts = 0;
-        forAll(pTimes, pi)
+        const labelList& cp = cellPoints[celli];
+        vector cellNi = vector::zero;
+        const point& cellCentre = cellCentres[celli];
+        forAll(cp, pointI)
         {
-            const label oldEdgeSign =
-                sign(pTimes[(pi + 1) % nPoints] - pTimes[pi]);
-            const label newEdgeSign =
-                sign(pTimes[(pi + 2) % nPoints] - pTimes[(pi + 1) % nPoints]);
-
-            if (newEdgeSign != oldEdgeSign)
-            {
-                nShifts++;
-            }
+            point vertex = points[cp[pointI]];
+            scalar w = 1.0/mag(vertex - cellCentre);
+            cellNi += w*vertexN[cp[pointI]];
         }
-
-        if (nShifts == 2)
-        {
-            dVf =
-                phi/magSf
-               *isoCutFace_.timeIntegratedArea(fPts, pTimes, dt, magSf, Un0);
-        }
-        else if (nShifts > 2)
-        {
-            // Triangle decompose the face
-            pointField fPts_tri(3);
-            scalarField pTimes_tri(3);
-            fPts_tri[0] = mesh_.faceCentres()[facei];
-            pTimes_tri[0] = ((fPts_tri[0] - x0) & n0)/Un0;
-            for (label pi = 0; pi < nPoints; pi++)
-            {
-                fPts_tri[1] = fPts[pi];
-                pTimes_tri[1] = pTimes[pi];
-                fPts_tri[2] = fPts[(pi + 1) % nPoints];
-                pTimes_tri[2] = pTimes[(pi + 1) % nPoints];
-                const scalar magSf_tri =
-                    mag
-                    (
-                        0.5
-                       *(fPts_tri[2] - fPts_tri[0])
-                       ^(fPts_tri[1] - fPts_tri[0])
-                    );
-                const scalar phi_tri = phi*magSf_tri/magSf;
-                dVf +=
-                    phi_tri
-                   /magSf_tri
-                   *isoCutFace_.timeIntegratedArea
-                    (
-                        fPts_tri,
-                        pTimes_tri,
-                        dt,
-                        magSf_tri,
-                        Un0
-                    );
-            }
-        }
-        else
-        {
-            if (debug)
-            {
-                WarningInFunction
-                    << "Warning: nShifts = " << nShifts << " on face " << facei
-                    << " with pTimes = " << pTimes << " owned by cell "
-                    << mesh_.faceOwner()[facei] << endl;
-            }
-        }
-
-        return dVf;
-    }
-    else
-    {
-        // Un0 is almost zero and isoFace is treated as stationary
-        isoCutFace_.calcSubFace(facei, f0);
-        const scalar alphaf = mag(isoCutFace_.subFaceArea()/magSf);
-
-        if (debug)
-        {
-            WarningInFunction
-                << "Un0 is almost zero (" << Un0
-                << ") - calculating dVf on face " << facei
-                << " using subFaceFraction giving alphaf = " << alphaf
-                << endl;
-        }
-
-        return phi*dt*alphaf;
+        cellNIn[celli] = cellNi/(mag(cellNi) + SMALL);
     }
 }
 
@@ -602,19 +478,20 @@ void Foam::isoAdvection::limitFluxes()
     // Get time step size
     const scalar dt = mesh_.time().deltaT().value();
 
-//    scalarField alphaNew = alpha1In_ - fvc::surfaceIntegrate(dVf_);
+    volScalarField alphaNew = alpha1_ - fvc::surfaceIntegrate(dVf_);
     const scalar aTol = 1.0e-12;          // Note: tolerances
-    const scalar maxAlphaMinus1 = 1;      // max(alphaNew - 1);
-    const scalar minAlpha = -1;           // min(alphaNew);
+    scalar maxAlphaMinus1 = gMax(alphaNew) - 1;      // max(alphaNew - 1);
+    scalar minAlpha = gMin(alphaNew);           // min(alphaNew);
     const label nUndershoots = 20;        // sum(neg0(alphaNew + aTol));
     const label nOvershoots = 20;         // sum(pos0(alphaNew - 1 - aTol));
     cellIsBounded_ = false;
 
+    Info << "isoAdvection: Before conservative bounding: min(alpha) = "
+        << minAlpha << ", max(alpha) = 1 + " << maxAlphaMinus1 << endl;
+
     // Loop number of bounding steps
     for (label n = 0; n < nAlphaBounds_; n++)
     {
-        Info<< "isoAdvection: bounding iteration " << n + 1 << endl;
-
         if (maxAlphaMinus1 > aTol) // Note: tolerances
         {
             DebugInfo << "Bound from above... " << endl;
@@ -974,9 +851,6 @@ void Foam::isoAdvection::advect()
     // Do the isoAdvection on surface cells
     timeIntegratedFlux();
 
-    // Synchronize processor patches
-    syncProcPatches(dVf_, phi_);
-
     // Adjust dVf for unbounded cells
     limitFluxes();
 
@@ -984,6 +858,11 @@ void Foam::isoAdvection::advect()
     alpha1_ -= fvc::surfaceIntegrate(dVf_);
     alpha1_.correctBoundaryConditions();
 
+    scalar maxAlphaMinus1 = gMax(alpha1In_) - 1;
+    scalar minAlpha = gMin(alpha1In_);
+    Info << "isoAdvection: After conservative bounding: min(alpha) = "
+        << minAlpha << ", max(alpha) = 1 + " << maxAlphaMinus1 << endl;
+
     // Apply non-conservative bounding mechanisms (clipping and snapping)
     // Note: We should be able to write out alpha before this is done!
     applyBruteForceBounding();
@@ -993,6 +872,9 @@ void Foam::isoAdvection::advect()
     writeBoundedCells();
 
     advectionTime_ += (mesh_.time().elapsedCpuTime() - advectionStartTime);
+    Info << "isoAdvection: time consumption = "
+        << label(100*advectionTime_/(mesh_.time().elapsedCpuTime() + SMALL))
+        << "%" << endl;
 }
 
 
@@ -1027,6 +909,8 @@ void Foam::isoAdvection::applyBruteForceBounding()
 
 void Foam::isoAdvection::writeSurfaceCells() const
 {
+    if (!mesh_.time().writeTime()) return;
+
     if (dict_.lookupOrDefault("writeSurfCells", false))
     {
         cellSet cSet
@@ -1049,6 +933,8 @@ void Foam::isoAdvection::writeSurfaceCells() const
 
 void Foam::isoAdvection::writeBoundedCells() const
 {
+    if (!mesh_.time().writeTime()) return;
+
     if (dict_.lookupOrDefault("writeBoundedCells", false))
     {
         cellSet cSet
@@ -1077,6 +963,9 @@ void Foam::isoAdvection::writeIsoFaces
     const DynamicList<List<point>>& faces
 ) const
 {
+
+    if (!writeIsoFacesToFile_ || !mesh_.time().writeTime()) return;
+
     // Writing isofaces to obj file for inspection, e.g. in paraview
     const fileName dirName
     (
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H
index 3e9db7f573..7aa04e6aec 100644
--- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H
@@ -194,18 +194,19 @@ class isoAdvection
             //- 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
+            //- Set ap_ values of celli's vertices in accordance with the
+            //  unit normal of celli as obtained from cellNoramlsIn.
+            void setCellVertexValues
             (
-                const label facei,
-                const vector& x0,
-                const vector& n0,
-                const scalar Un0,
-                const scalar f0,
-                const scalar dt,
-                const scalar phi,
-                const scalar magSf
+                const label celli,
+                const vectorField& cellNormalsIn
+            );
+
+            //- Function used to normalise and smoothen grad(alpha) in case
+            //  gradAlphaBasedNormal_ is true.
+            void normaliseAndSmooth
+            (
+                volVectorField& cellN
             );
 
             //- For a given cell return labels of faces fluxing out of this cell
@@ -255,6 +256,17 @@ class isoAdvection
                 bsn0_.clear();
                 bsUn0_.clear();
                 bsf0_.clear();
+
+                if (mesh_.topoChanging())
+                {
+                    // Introduced resizing to cope with changing meshes
+                    checkBounding_.resize(mesh_.nCells());
+                    cellIsBounded_.resize(mesh_.nCells());
+                    ap_.resize(mesh_.nPoints());
+                }
+                checkBounding_ = false;
+                cellIsBounded_ = false;
+
             }
 
         // Face value functions needed for random face access where the face
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C
index 30aaa5f72b..8e26f5c606 100644
--- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C
@@ -327,7 +327,7 @@ Foam::label Foam::isoCutCell::calcSubCell
     const scalar isoValue
 )
 {
-    // Populate isoCutFaces_, isoCutFacePoints_, fullySubFaces_, isoFaceCentres_
+    // Populate isoCutFaces_, isoCutFacePoints_, fullySubFaces_, isoFaceCentre_
     // and isoFaceArea_.
 
     clearStorage();
@@ -501,7 +501,7 @@ Foam::label Foam::isoCutCell::vofCutCell
         << "vofCutCell for cell " << celli << " with alpha1 = "
         << alpha1 << " ------" << endl;
 
-    // Finding cell vertex extrema values
+    // Finding cell vertex extremum values
     const labelList& pLabels = mesh_.cellPoints(celli);
     scalarField fvert(pLabels.size());
     forAll(pLabels, pi)
@@ -574,15 +574,17 @@ Foam::label Foam::isoCutCell::vofCutCell
     calcSubCell(celli, f3);
     a3 = volumeOfFluid();
 
-    scalar f4 = f1 + 2*(f2 - f1)/3;
+    scalar f4 = f1 + (f2 - f1)*scalar(2)/scalar(3);
     calcSubCell(celli, f4);
     scalar a4 = volumeOfFluid();
 
     // Building and solving Vandermonde matrix equation
     scalarField a(4), f(4), C(4);
     {
-        a[0] = a1, a[1] = a3, a[2] = a4, a[3] = a2;
-        f[0] = 0, f[1] = (f3-f1)/(f2-f1), f[2] = (f4-f1)/(f2-f1), f[3] = 1;
+        a[0] = a1, f[0] = 0;
+        a[1] = a3, f[1] = (f3 - f1)/(f2 - f1);
+        a[2] = a4, f[2] = (f4 - f1)/(f2 - f1);
+        a[3] = a2, f[3] = 1;
         scalarSquareMatrix M(4);
         forAll(f, i)
         {
@@ -598,7 +600,6 @@ Foam::label Foam::isoCutCell::vofCutCell
     }
 
     // Finding root with Newton method
-
     f3 = f[1]; a3 = a[1];
     label nIter = 0;
     scalar res = mag(a3 - alpha1);
@@ -639,7 +640,7 @@ Foam::label Foam::isoCutCell::vofCutCell
 
     // If tolerance not met use the secant method  with f3 as a hopefully very
     // good initial guess to crank res the last piece down below tol
-
+    // Note: This is expensive because subcell is recalculated every iteration
     scalar x2 = f3;
     scalar g2 = VOF - alpha1;
     scalar x1 = max(1e-3*(f2 - f1), 100*SMALL);
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
index 83292012a0..43bc4bb610 100644
--- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
@@ -338,6 +338,152 @@ void Foam::isoCutFace::clearStorage()
 }
 
 
+Foam::scalar Foam::isoCutFace::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
+)
+{
+/* Temporarily taken out
+    // Treating rare cases where isoface normal is not calculated properly
+    if (mag(n0) < 0.5)
+    {
+        scalar alphaf = 0;
+        scalar waterInUpwindCell = 0;
+
+        if (phi > 10*SMALL || !mesh_.isInternalFace(facei))
+        {
+            const label upwindCell = mesh_.faceOwner()[facei];
+            alphaf = alpha1In_[upwindCell];
+            waterInUpwindCell = alphaf*mesh_.cellVolumes()[upwindCell];
+        }
+        else
+        {
+            const label upwindCell = mesh_.faceNeighbour()[facei];
+            alphaf = alpha1In_[upwindCell];
+            waterInUpwindCell = alphaf*mesh_.cellVolumes()[upwindCell];
+        }
+
+        if (debug)
+        {
+            WarningInFunction
+                << "mag(n0) = " << mag(n0)
+                << " so timeIntegratedFlux calculates dVf from upwind"
+                << " cell alpha value: " << alphaf << endl;
+        }
+
+        return min(alphaf*phi*dt, waterInUpwindCell);
+    }
+*/
+
+    // Find sorted list of times where the isoFace will arrive at face points
+    // given initial position x0 and velocity Un0*n0
+
+    // Get points for this face
+    const face& f = mesh_.faces()[facei];
+    const pointField fPts(f.points(mesh_.points()));
+    const label nPoints = fPts.size();
+
+    scalarField pTimes(fPts.size());
+    if (mag(Un0) > 10*SMALL) // Note: tolerances
+    {
+        // Here we estimate time of arrival to the face points from their normal
+        // distance to the initial surface and the surface normal velocity
+
+        pTimes = ((fPts - x0) & n0)/Un0;
+
+        scalar dVf = 0;
+
+        // Check if pTimes changes direction more than twice when looping face
+        label nShifts = 0;
+        forAll(pTimes, pi)
+        {
+            const label oldEdgeSign =
+                sign(pTimes[(pi + 1) % nPoints] - pTimes[pi]);
+            const label newEdgeSign =
+                sign(pTimes[(pi + 2) % nPoints] - pTimes[(pi + 1) % nPoints]);
+
+            if (newEdgeSign != oldEdgeSign)
+            {
+                nShifts++;
+            }
+        }
+
+        if (nShifts == 2)
+        {
+            dVf = phi/magSf*timeIntegratedArea(fPts, pTimes, dt, magSf, Un0);
+        }
+        else if (nShifts > 2)
+        {
+            // Triangle decompose the face
+            pointField fPts_tri(3);
+            scalarField pTimes_tri(3);
+            fPts_tri[0] = mesh_.faceCentres()[facei];
+            pTimes_tri[0] = ((fPts_tri[0] - x0) & n0)/Un0;
+            for (label pi = 0; pi < nPoints; pi++)
+            {
+                fPts_tri[1] = fPts[pi];
+                pTimes_tri[1] = pTimes[pi];
+                fPts_tri[2] = fPts[(pi + 1) % nPoints];
+                pTimes_tri[2] = pTimes[(pi + 1) % nPoints];
+                const scalar magSf_tri =
+                    mag
+                    (
+                        0.5
+                       *(fPts_tri[2] - fPts_tri[0])
+                       ^(fPts_tri[1] - fPts_tri[0])
+                    );
+                const scalar phi_tri = phi*magSf_tri/magSf;
+                dVf += phi_tri/magSf_tri
+                   *timeIntegratedArea
+                    (
+                        fPts_tri,
+                        pTimes_tri,
+                        dt,
+                        magSf_tri,
+                        Un0
+                    );
+            }
+        }
+        else
+        {
+            if (debug)
+            {
+                WarningInFunction
+                    << "Warning: nShifts = " << nShifts << " on face " << facei
+                    << " with pTimes = " << pTimes << " owned by cell "
+                    << mesh_.faceOwner()[facei] << endl;
+            }
+        }
+
+        return dVf;
+    }
+    else
+    {
+        // Un0 is almost zero and isoFace is treated as stationary
+        calcSubFace(facei, f0);
+        const scalar alphaf = mag(subFaceArea()/magSf);
+
+        if (debug)
+        {
+            WarningInFunction
+                << "Un0 is almost zero (" << Un0
+                << ") - calculating dVf on face " << facei
+                << " using subFaceFraction giving alphaf = " << alphaf
+                << endl;
+        }
+
+        return phi*dt*alphaf;
+    }
+}
+
+
 Foam::scalar Foam::isoCutFace::timeIntegratedArea
 (
     const pointField& fPts,
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.H b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.H
index 7afbbf3639..46e746977f 100644
--- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.H
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.H
@@ -159,7 +159,21 @@ public:
         const DynamicList<point>& surfacePoints() const;
 
         void clearStorage();
-
+        
+        //- 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
+        );
+        
         //- Calculate time integrated area for a face
         scalar timeIntegratedArea
         (
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/U b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/U
new file mode 100644
index 0000000000..5ac9838bc8
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/U
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    atmosphere
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+    walls
+    {
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/alpha.water b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/alpha.water
new file mode 100644
index 0000000000..33c24db83b
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/alpha.water
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    walls
+    {
+        type            zeroGradient;
+    }
+
+    obstacle
+    {
+        type            zeroGradient;
+    }
+
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/p_rgh b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/p_rgh
new file mode 100644
index 0000000000..481f0a268f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/0.orig/p_rgh
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    walls
+    {
+        type            fixedFluxPressure;
+        phi             phiAbs;
+        value           uniform 0;
+    }
+
+    obstacle
+    {
+        type            fixedFluxPressure;
+        phi             phiAbs;
+        value           uniform 0;
+    }
+
+    atmosphere
+    {
+        type            totalPressure;
+        p0              uniform 0;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/Allclean b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/Allclean
new file mode 100755
index 0000000000..d1295e344e
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/Allclean
@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+rm system/cellSetDict > /dev/null 2>&1
+rm -rf 0 > /dev/null 2>&1
+
+cleanCase
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/Allrun b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/Allrun
new file mode 100755
index 0000000000..f59721c39c
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/Allrun
@@ -0,0 +1,17 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+cp -r 0.orig 0 > /dev/null 2>&1
+runApplication blockMesh
+#runApplication setSet -batch createObstacle.setSet
+runApplication topoSet
+runApplication subsetMesh -overwrite c0 -patch walls
+runApplication setFields
+#runApplication decomposePar
+#runParallel `getApplication`
+runApplication `getApplication`
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/dynamicMeshDict b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/dynamicMeshDict
new file mode 100644
index 0000000000..7acf1ce97a
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/dynamicMeshDict
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicRefineFvMesh;
+
+// How often to refine
+refineInterval  1;
+
+// Field to be refinement on
+field           alpha.water;
+
+// Refine field inbetween lower..upper
+lowerRefineLevel 0.001;
+upperRefineLevel 0.999;
+
+// If value < unrefineLevel unrefine
+unrefineLevel   10;
+
+// Have slower than 2:1 refinement
+nBufferLayers   1;
+
+// Refine cells only up to maxRefinement levels
+maxRefinement   2;
+
+// Stop refinement if maxCells reached
+maxCells        15000000;
+
+// Flux field and corresponding velocity field. Fluxes on changed
+// faces get recalculated by interpolating the velocity. Use 'none'
+// on surfaceScalarFields that do not need to be reinterpolated.
+correctFluxes
+(
+    (phi none)
+    (nHatf none)
+    (rhoPhi none)
+    (alphaPhi none)
+    (ghf none)
+    (phi0 none)
+    (dVf_ none)
+);
+
+// Write the refinement level as a volScalarField
+dumpLevel       true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/g b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/g
new file mode 100644
index 0000000000..ac04af46dd
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 0 -9.81);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/transportProperties b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/transportProperties
new file mode 100644
index 0000000000..37551da4d5
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/transportProperties
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+water
+{
+    transportModel  Newtonian;
+    nu              1e-06;
+    rho             1000;
+}
+
+air
+{
+    transportModel  Newtonian;
+    nu              1.48e-05;
+    rho             1;
+}
+
+sigma           0.07;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/turbulenceProperties b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/turbulenceProperties
new file mode 100644
index 0000000000..8b144cad5f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/blockMeshDict b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/blockMeshDict
new file mode 100644
index 0000000000..7cefb52001
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/blockMeshDict
@@ -0,0 +1,68 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 1)
+    (1 0 1)
+    (1 1 1)
+    (0 1 1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (16 16 16) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    atmosphere
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    walls
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+            (2 6 5 1)
+            (1 5 4 0)
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/controlDict b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/controlDict
new file mode 100644
index 0000000000..5e9079b14e
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/controlDict
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     interIsoFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime; //timeStep;
+
+writeInterval   0.02;//1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.75;
+maxAlphaCo      0.75;
+maxDeltaT       1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/decomposeParDict b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/decomposeParDict
new file mode 100644
index 0000000000..b2a4b7e2df
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+method              simple;
+
+simpleCoeffs
+{
+    n               (2 2 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (2 2 1);
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/fvSchemes b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/fvSchemes
new file mode 100644
index 0000000000..f808a3220e
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/fvSchemes
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(rhoPhi,U)   Gauss upwind;
+    div(phi,alpha)  Gauss vanLeer;
+    div(phirb,alpha) Gauss linear;
+    div(phi,k)      Gauss upwind;
+    div(phi,omega)  Gauss upwind;
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+wallDist
+{
+    method meshWave;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/fvSolution b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/fvSolution
new file mode 100644
index 0000000000..fda530bd5b
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/fvSolution
@@ -0,0 +1,89 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha.water.*"
+    {
+        isoFaceTol      1e-6;
+        surfCellTol     1e-6;
+        nAlphaBounds    3;
+        snapTol         1e-12;
+        clip            true;
+        writeSurfCells  false;
+        writeBoundedCells false;
+        writeIsoFaces   false;
+
+        nAlphaCorr      1;
+        nAlphaSubCycles 1;
+        cAlpha          1;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.01;
+        smoother        DIC;
+        cacheAgglomeration no;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+        tolerance       1e-9;
+    }
+
+    "pcorr.*"
+    {
+        $p_rghFinal;
+//        tolerance       0.0001;
+        tolerance       1e-08;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-07;
+        relTol          0;
+        nSweeps         1;
+    }
+
+    "(k|omega|B|nuTilda).*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-08;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor no;
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 1;
+    correctPhi     yes;
+
+    pRefPoint      (0.51 0.51 0.51);
+    pRefValue      0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/setFieldsDict b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/setFieldsDict
new file mode 100644
index 0000000000..7a8e532cce
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/setFieldsDict
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.water 0
+    volVectorFieldValue U (0 0 0)
+);
+
+regions
+(
+    boxToCell
+    {
+        box (0 0 0) (0.6 0.1875 0.75);
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 1
+        );
+    }
+ );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/topoSetDict b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/topoSetDict
new file mode 100644
index 0000000000..d23d2fc637
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/topoSetDict
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    c0;
+        type    cellSet;
+        action  clear;
+    }
+    {
+        name    c0;
+        type    cellSet;
+        action  invert;
+    }
+    {
+        name    c0;
+        type    cellSet;
+        action  delete;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (0.375 0.375 0) (0.625 0.625 0.25);
+        }
+    }
+);
+
+// ************************************************************************* //