Merge branch 'feature-isoAdvector-AMR' into 'develop-pre-release'

Feature iso advector AMR

See merge request OpenFOAM-plus!205

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 \

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

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();
+    if (pimple.nCorrPIMPLE() > 1)
+    {
+        // If nOuterCorrectors > 1 then for all but the first loop the advection
+        // of alpha is done using an average, 0.5*phi+0.5*phiNew where phi is 
+        // the flux at the beginning of the time step and phiNew is the flux 
+        // estimate at the end of the time step from the previous outer 
+        // iteration. Similarly we use 0.5*U + 0.5*UNew in later iterations.
+        if (pimple.firstIter())
+        {
+            // To recalculate the alpha1 update in subsequent iterations, we 
+            // must store its current value before overwriting with the new 
+            // value
+            alpha1.prevIter();
+            // Storing initial phi and U for use in later outer iterations.
+            phi.storePrevIter();
+            U.storePrevIter();
+        }
+        else
+        {
+            // Resetting alpha1 to value before advection in first PIMPLE 
+            // iteration.
+            alpha1 = alpha1.prevIter();
+            // Setting U and phi with which to advect interface.
+            U = 0.5*U.prevIter() + 0.5*U;
+            phi = 0.5*phi.prevIter() + 0.5*phi;
+        }
+    }
+    
+    // Updating alpha1
+    advector.advect();
+    #include "rhofs.H"
+    rhoPhi = advector.getRhoPhi(rho1f, rho2f);
 
-    // 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();
+    if (!pimple.firstIter())
+    {
+        // Resetting U and phi to value at latest iteration.
+        U = 2.0*U - U.prevIter();
+        phi = 2.0*phi - phi.prevIter();
+    }
 
-// Update rhoPhi
-rhoPhi = advector.getRhoPhi(rho1, rho2);
+    alpha2 = 1.0 - alpha1;
+    mixture.correct();
 
-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
+    << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
     << endl;

applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H

@@ -13,8 +13,6 @@ if (nAlphaSubCycles > 1)
         dimensionedScalar(rhoPhi.dimensions(), Zero)
     );
 
-    tmp<volScalarField> trSubDeltaT;
-
     for
     (
         subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);

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
 );

applications/solvers/multiphase/interIsoFoam/createFields.H

+#include "createRDeltaT.H"
+
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
 (
@@ -120,4 +122,6 @@ mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
 #include "createMRF.H"
-#include "createIsoAdvection.H"
+#include "createFvOptions.H"
+
+isoAdvection advector(alpha1, phi, U);

applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H

-// Defining isoAdvection
-isoAdvection advector(alpha1, phi, U);

applications/solvers/multiphase/interIsoFoam/initCorrectPhi.H

+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"
+}

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,28 +71,24 @@ 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 "initCorrectPhi.H"
+    #include "createUfIfPresent.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 "readDyMControls.H"
         #include "CourantNo.H"
         #include "alphaCourantNo.H"
         #include "setDeltaT.H"
@@ -104,6 +100,39 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                mesh.update();
+
+                if (mesh.changing())
+                {
+
+                    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"
 

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));
+    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();
+    }
 }

applications/solvers/multiphase/interIsoFoam/rhofs.H

+const dimensionedScalar& rho1f(rho1);
+const dimensionedScalar& rho2f(rho2);

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

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

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                 isoAdvector | Copyright (C) 2016-2017 DHI
+              Modified work | Copyright (C) 2018 Johan Roenby
 -------------------------------------------------------------------------------
 
 License
@@ -194,18 +195,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 +257,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

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                 isoAdvector | Copyright (C) 2016-2017 DHI
+              Modified work | Copyright (C) 2018 Johan Roenby
 -------------------------------------------------------------------------------
 
 License
@@ -327,7 +328,7 @@ Foam::label Foam::isoCutCell::calcSubCell
     const scalar isoValue
 )
 {
-    // Populate isoCutFaces_, isoCutFacePoints_, fullySubFaces_, isoFaceCentres_
+    // Populate isoCutFaces_, isoCutFacePoints_, fullySubFaces_, isoFaceCentre_
     // and isoFaceArea_.
 
     clearStorage();
@@ -501,7 +502,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 +575,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 +601,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 +641,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);

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

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                 isoAdvector | Copyright (C) 2016-2017 DHI
+              Modified work | Copyright (C) 2018 Johan Roenby
 -------------------------------------------------------------------------------
 
 License

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                 isoAdvector | Copyright (C) 2016-2017 DHI
+              Modified work | Copyright (C) 2018 Johan Roenby
 -------------------------------------------------------------------------------
 
 License
@@ -338,6 +339,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;