functionObjects::scalarTransport: Added support for MULES with sub-cycling and semi-implicit options

Description
    Evolves a passive scalar transport equation.

    - To specify the field name set the \c field entry
    - To employ the same numerical schemes as another field set
      the \c schemesField entry,
    - The \c diffusivity entry can be set to \c none, \c constant, \c viscosity
    - A constant diffusivity is specified with the \c D entry,
    - If a momentum transport model is available and the \c viscosity
      diffusivety option specified an effective diffusivity may be constructed
      from the laminar and turbulent viscosities using the diffusivity
      coefficients \c alphal and \c alphat:
      \verbatim
          D = alphal*nu + alphat*nut
      \endverbatim

    Example:
    \verbatim
        #includeFunc scalarTransport(T, alphaD=1, alphaDt=1)
    \endverbatim

    For incompressible flow the passive scalar may optionally be solved with the
    MULES limiter and sub-cycling or semi-implicit in order to maintain
    boundedness, particularly if a compressive, PLIC or MPLIC convection
    scheme is used.

    Example:
    \verbatim
        #includeFunc scalarTransport(tracer, diffusion=none)

    with scheme specification:
        div(phi,tracer)     Gauss interfaceCompression vanLeer 1;

    and solver specification:
        tracer
        {
            nCorr      1;
            nSubCycles 3;

            MULESCorr       no;
            nLimiterIter    5;
            applyPrevCorr   yes;

            solver          smoothSolver;
            smoother        symGaussSeidel;
            tolerance       1e-8;
            relTol          0;

            diffusion
            {
                solver          smoothSolver;
                smoother        symGaussSeidel;
                tolerance       1e-8;
                relTol          0;
            }
        }
    \endverbatim

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -41,6 +41,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "interfaceCompression.H"
 #include "CMULES.H"
 #include "EulerDdtScheme.H"
 #include "localEulerDdtScheme.H"

applications/solvers/multiphase/interFoam/interFoam.C

@@ -33,6 +33,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "interfaceCompression.H"
 #include "CMULES.H"
 #include "EulerDdtScheme.H"
 #include "localEulerDdtScheme.H"

applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C

@@ -32,6 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "interfaceCompression.H"
 #include "CMULES.H"
 #include "localEulerDdtScheme.H"
 #include "subCycle.H"

etc/caseDicts/postProcessing/solvers/scalarTransport

@@ -21,5 +21,9 @@ Description
 field           <fieldName>;    // Name of the transported scalar
 schemesField    $field;         // Name of the field from which to use schemes
                                 // and solvers settings
+diffusion       viscosity;
+
+alphal          1;
+alphat          1;
 
 // ************************************************************************* //

src/functionObjects/solvers/Make/options

@@ -4,6 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
@@ -13,6 +14,7 @@ LIB_LIBS = \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
     -lcompressibleMomentumTransportModels \
+    -ltwoPhaseMixture \
     -lspecie \
     -lfiniteVolume \
     -lmeshTools \

src/functionObjects/solvers/scalarTransport/scalarTransport.C

@@ -26,13 +26,23 @@ License
 #include "scalarTransport.H"
 #include "surfaceFields.H"
 #include "fvmDdt.H"
+#include "fvcDdt.H"
 #include "fvmDiv.H"
 #include "fvmLaplacian.H"
 #include "fvmSup.H"
+#include "fvcFlux.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 #include "incompressibleMomentumTransportModel.H"
 #include "compressibleMomentumTransportModel.H"
+
+#include "CMULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
+#include "subCycle.H"
+#include "interfaceCompression.H"
+
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -53,16 +63,33 @@ namespace functionObjects
 }
 
 
+template<>
+const char* Foam::NamedEnum
+<
+    Foam::functionObjects::scalarTransport::diffusionType,
+    3
+>::names[] =
+{
+    "none",
+    "constant",
+    "viscosity"
+};
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::scalarTransport::diffusionType,
+    3
+> Foam::functionObjects::scalarTransport::diffusionTypeNames_;
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
+Foam::tmp<Foam::volScalarField>
-(
+Foam::functionObjects::scalarTransport::D() const
-    const surfaceScalarField& phi
-) const
 {
     const word Dname("D" + s_.name());
 
-    if (constantD_)
+    if (diffusion_ == diffusionType::constant)
     {
         return volScalarField::New
         (
@@ -71,13 +98,7 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
             dimensionedScalar(Dname, dimViscosity, D_)
         );
     }
-    else if
+    else
-    (
-        mesh_.foundObject<momentumTransportModel>
-        (
-            momentumTransportModel::typeName
-        )
-    )
     {
         const momentumTransportModel& turbulence =
             mesh_.lookupObject<momentumTransportModel>
@@ -85,16 +106,7 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
                 momentumTransportModel::typeName
             );
 
-        return alphaD_*turbulence.nu() + alphaDt_*turbulence.nut();
+        return alphal_*turbulence.nu() + alphat_*turbulence.nut();
-    }
-    else
-    {
-        return volScalarField::New
-        (
-            Dname,
-            mesh_,
-            dimensionedScalar(Dname, phi.dimensions()/dimLength, 0)
-        );
     }
 }
 
@@ -110,6 +122,7 @@ Foam::functionObjects::scalarTransport::scalarTransport
 :
     fvMeshFunctionObject(name, runTime, dict),
     fieldName_(dict.lookupOrDefault<word>("field", "s")),
+    diffusion_(diffusionType::none),
     D_(0),
     nCorr_(0),
     s_
@@ -123,9 +136,49 @@ Foam::functionObjects::scalarTransport::scalarTransport
             IOobject::AUTO_WRITE
         ),
         mesh_
-    )
+    ),
+    MULES_(false),
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mesh_.V()), 1.0/3.0)
+    ),
+    sRestart_(false)
 {
     read(dict);
+
+    const dictionary& controls = mesh_.solverDict(s_.name());
+
+    if (controls.found("nSubCycles"))
+    {
+        MULES_ = true;
+
+        typeIOobject<surfaceScalarField> sPhiHeader
+        (
+            IOobject::groupName("sPhi", s_.group()),
+            runTime.timeName(),
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        );
+
+        sRestart_ = sPhiHeader.headerOk();
+
+        if (sRestart_)
+        {
+            Info << "Restarting s" << endl;
+        }
+
+        const surfaceScalarField& phi =
+            mesh_.lookupObject<surfaceScalarField>(phiName_);
+
+        // Scalar volumetric flux
+        tsPhi_ = new surfaceScalarField
+        (
+            sPhiHeader,
+            phi*fvc::interpolate(s_)
+        );
+    }
 }
 
 
@@ -145,9 +198,22 @@ bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
     rhoName_ = dict.lookupOrDefault<word>("rho", "rho");
     schemesField_ = dict.lookupOrDefault<word>("schemesField", fieldName_);
 
-    constantD_ = dict.readIfPresent("D", D_);
+    diffusion_ = diffusionTypeNames_.read(dict.lookup("diffusion"));
-    alphaD_ = dict.lookupOrDefault("alphaD", 1.0);
+
-    alphaDt_ = dict.lookupOrDefault("alphaDt", 1.0);
+    switch(diffusion_)
+    {
+        case diffusionType::none:
+            break;
+
+        case diffusionType::constant:
+            dict.lookup("D") >> D_;
+            break;
+
+        case diffusionType::viscosity:
+            dict.lookup("alphal") >> alphal_;
+            dict.lookup("alphat") >> alphat_;
+            break;
+    }
 
     dict.readIfPresent("nCorr", nCorr_);
 
@@ -168,11 +234,7 @@ bool Foam::functionObjects::scalarTransport::execute()
     const surfaceScalarField& phi =
         mesh_.lookupObject<surfaceScalarField>(phiName_);
 
-    // Calculate the diffusivity
+    const word divScheme("div(phi," + schemesField_ + ")");
-    volScalarField D("D" + s_.name(), this->D(phi));
-
-    word divScheme("div(phi," + schemesField_ + ")");
-    word laplacianScheme("laplacian(" + D.name() + "," + schemesField_ + ")");
 
     // Set under-relaxation coeff
     scalar relaxCoeff = 0.0;
@@ -189,25 +251,38 @@ bool Foam::functionObjects::scalarTransport::execute()
 
     if (phi.dimensions() == dimVolume/dimTime)
     {
-        for (int i=0; i<=nCorr_; i++)
+        if (MULES_)
         {
-            fvScalarMatrix sEqn
+            subCycleMULES();
-            (
-                fvm::ddt(s_)
-              + fvm::div(phi, s_, divScheme)
-              - fvm::laplacian(D, s_, laplacianScheme)
-             ==
-                fvModels.source(s_)
-            );
-
-            sEqn.relax(relaxCoeff);
-
-            fvConstraints.constrain(sEqn);
-
-            sEqn.solve(schemesField_);
 
             fvConstraints.constrain(s_);
         }
+        else
+        {
+            for (int i=0; i<=nCorr_; i++)
+            {
+                fvScalarMatrix sEqn
+                (
+                    fvm::ddt(s_)
+                  + fvm::div(phi, s_, divScheme)
+                 ==
+                    fvModels.source(s_)
+                );
+
+                if (diffusion_ != diffusionType::none)
+                {
+                    sEqn -= fvm::laplacian(D(), s_);
+                }
+
+                sEqn.relax(relaxCoeff);
+
+                fvConstraints.constrain(sEqn);
+
+                sEqn.solve(schemesField_);
+
+                fvConstraints.constrain(s_);
+            }
+        }
     }
     else if (phi.dimensions() == dimMass/dimTime)
     {
@@ -220,11 +295,15 @@ bool Foam::functionObjects::scalarTransport::execute()
             (
                 fvm::ddt(rho, s_)
               + fvm::div(phi, s_, divScheme)
-              - fvm::laplacian(rho*D, s_, laplacianScheme)
              ==
                 fvModels.source(rho, s_)
             );
 
+            if (diffusion_ != diffusionType::none)
+            {
+                sEqn -= fvm::laplacian(rho*D(), s_);
+            }
+
             sEqn.relax(relaxCoeff);
 
             fvConstraints.constrain(sEqn);
@@ -248,6 +327,303 @@ bool Foam::functionObjects::scalarTransport::execute()
 }
 
 
+void Foam::functionObjects::scalarTransport::subCycleMULES()
+{
+    const dictionary& controls = mesh_.solverDict(s_.name());
+    const label nSubCycles(controls.lookup<label>("nSubCycles"));
+    const bool LTS = fv::localEulerDdt::enabled(mesh_);
+
+    if (nSubCycles > 1)
+    {
+        tmp<volScalarField> trSubDeltaT;
+
+        if (LTS)
+        {
+            trSubDeltaT =
+                fv::localEulerDdt::localRSubDeltaT(mesh_, nSubCycles);
+        }
+
+        for
+        (
+            subCycle<volScalarField> sSubCycle(s_, nSubCycles);
+            !(++sSubCycle).end();
+        )
+        {
+            solveMULES();
+        }
+    }
+    else
+    {
+        solveMULES();
+    }
+
+
+    // Apply the diffusion term separately to allow implicit solution
+    // and boundedness of the explicit advection
+    if (diffusion_ != diffusionType::none)
+    {
+        fvScalarMatrix sEqn
+        (
+            fvm::ddt(s_) - fvc::ddt(s_)
+          - fvm::laplacian(D(), s_)
+        );
+
+        sEqn.solve(controls.subDict("diffusion"));
+
+        Info<< s_.name() << " volume fraction = "
+            << s_.weightedAverage(mesh_.V()).value()
+            << "  Min(" << s_.name() << ") = " << min(s_).value()
+            << "  Max(" << s_.name() << ") = " << max(s_).value()
+            << endl;
+    }
+}
+
+
+void Foam::functionObjects::scalarTransport::solveMULES()
+{
+    const dictionary& controls = mesh_.solverDict(s_.name());
+    const label nCorr(controls.lookup<label>("nCorr"));
+    const label nSubCycles(controls.lookup<label>("nSubCycles"));
+    const bool MULESCorr(controls.lookupOrDefault<Switch>("MULESCorr", false));
+
+    // Apply the compression correction from the previous iteration
+    // Improves efficiency for steady-simulations but can only be applied
+    // once the s field is reasonably steady, i.e. fully developed
+    const bool applyPrevCorr
+    (
+        controls.lookupOrDefault<Switch>("applyPrevCorr", false)
+    );
+
+    const bool LTS = fv::localEulerDdt::enabled(mesh_);
+
+    const word divScheme("div(phi," + schemesField_ + ")");
+
+    const surfaceScalarField& phi =
+        mesh_.lookupObject<surfaceScalarField>(phiName_);
+
+    surfaceScalarField& sPhi_ = tsPhi_.ref();
+
+    const word sScheme(mesh_.divScheme(divScheme)[1].wordToken());
+
+    // If a compressive convection scheme is used
+    // the interface normal must be cached
+    tmp<surfaceScalarField> nHatf;
+
+    if (compressionSchemes.found(sScheme))
+    {
+        const surfaceVectorField gradsf(fvc::interpolate(fvc::grad(s_)));
+
+        nHatf = new surfaceScalarField
+        (
+            IOobject
+            (
+                "nHatf",
+                s_.time().timeName(),
+                mesh_
+            ),
+            gradsf/(mag(gradsf) + deltaN_) & mesh_.Sf()
+        );
+    }
+
+    // Set the off-centering coefficient according to ddt scheme
+    scalar ocCoeff = 0;
+    {
+        tmp<fv::ddtScheme<scalar>> tddtS
+        (
+            fv::ddtScheme<scalar>::New
+            (
+                mesh_,
+                mesh_.ddtScheme("ddt(s)")
+            )
+        );
+        const fv::ddtScheme<scalar>& ddtS = tddtS();
+
+        if
+        (
+            isType<fv::EulerDdtScheme<scalar>>(ddtS)
+         || isType<fv::localEulerDdtScheme<scalar>>(ddtS)
+        )
+        {
+            ocCoeff = 0;
+        }
+        else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtS))
+        {
+            if (nSubCycles > 1)
+            {
+                FatalErrorInFunction
+                    << "Sub-cycling is not supported "
+                       "with the CrankNicolson ddt scheme"
+                    << exit(FatalError);
+            }
+
+            if
+            (
+                sRestart_
+             || mesh_.time().timeIndex() > mesh_.time().startTimeIndex() + 1
+            )
+            {
+                ocCoeff =
+                    refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtS)
+                   .ocCoeff();
+            }
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Only Euler and CrankNicolson ddt schemes are supported"
+                << exit(FatalError);
+        }
+    }
+
+    // Set the time blending factor, 1 for Euler
+    scalar cnCoeff = 1.0/(1.0 + ocCoeff);
+
+    tmp<surfaceScalarField> phiCN(phi);
+
+    // Calculate the Crank-Nicolson off-centred volumetric flux
+    if (ocCoeff > 0)
+    {
+        phiCN = surfaceScalarField::New
+        (
+            "phiCN",
+            cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime()
+        );
+    }
+
+    if (MULESCorr)
+    {
+        fvScalarMatrix sEqn
+        (
+            (
+                LTS
+              ? fv::localEulerDdtScheme<scalar>(mesh_).fvmDdt(s_)
+              : fv::EulerDdtScheme<scalar>(mesh_).fvmDdt(s_)
+            )
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh_,
+                phiCN,
+                upwind<scalar>(mesh_, phiCN)
+            ).fvmDiv(phiCN, s_)
+        );
+
+        sEqn.solve();
+
+        Info<< s_.name() << " volume fraction = "
+            << s_.weightedAverage(mesh_.Vsc()).value()
+            << "  Min(" << s_.name() << ") = " << min(s_).value()
+            << "  Max(" << s_.name() << ") = " << max(s_).value()
+            << endl;
+
+        tmp<surfaceScalarField> tsPhiUD(sEqn.flux());
+        sPhi_ = tsPhiUD();
+
+        if (applyPrevCorr && tsPhiCorr0_.valid())
+        {
+            Info<< "Applying the previous iteration compression flux" << endl;
+            MULES::correct
+            (
+                geometricOneField(),
+                s_,
+                sPhi_,
+                tsPhiCorr0_.ref(),
+                oneField(),
+                zeroField()
+            );
+
+            sPhi_ += tsPhiCorr0_();
+        }
+
+        // Cache the upwind-flux
+        tsPhiCorr0_ = tsPhiUD;
+    }
+
+    for (int sCorr=0; sCorr<nCorr; sCorr++)
+    {
+        // Split operator
+        tmp<surfaceScalarField> tsPhiUn
+        (
+            fvc::flux
+            (
+                phiCN(),
+                (cnCoeff*s_ + (1.0 - cnCoeff)*s_.oldTime())(),
+                mesh_.divScheme(divScheme)
+            )
+        );
+
+        if (MULESCorr)
+        {
+            tmp<surfaceScalarField> tsPhiCorr(tsPhiUn() - sPhi_);
+            volScalarField s0("s0", s_);
+
+            MULES::correct
+            (
+                geometricOneField(),
+                s_,
+                tsPhiUn(),
+                tsPhiCorr.ref(),
+                oneField(),
+                zeroField()
+            );
+
+            // Under-relax the correction for all but the 1st corrector
+            if (sCorr == 0)
+            {
+                sPhi_ += tsPhiCorr();
+            }
+            else
+            {
+                s_ = 0.5*s_ + 0.5*s0;
+                sPhi_ += 0.5*tsPhiCorr();
+            }
+        }
+        else
+        {
+            sPhi_ = tsPhiUn;
+
+            MULES::explicitSolve
+            (
+                geometricOneField(),
+                s_,
+                phiCN,
+                sPhi_,
+                oneField(),
+                zeroField()
+            );
+        }
+    }
+
+    if (applyPrevCorr && MULESCorr)
+    {
+        tsPhiCorr0_ = sPhi_ - tsPhiCorr0_;
+        tsPhiCorr0_.ref().rename("sPhiCorr0");
+    }
+    else
+    {
+        tsPhiCorr0_.clear();
+    }
+
+    if
+    (
+        word(mesh_.ddtScheme("ddt(s)"))
+     == fv::CrankNicolsonDdtScheme<scalar>::typeName
+    )
+    {
+        if (ocCoeff > 0)
+        {
+            // Calculate the end-of-time-step s flux
+            sPhi_ = (sPhi_ - (1.0 - cnCoeff)*sPhi_.oldTime())/cnCoeff;
+        }
+    }
+
+    Info<< s_.name() << "volume fraction = "
+        << s_.weightedAverage(mesh_.Vsc()).value()
+        << "  Min(" << s_.name() << ") = " << min(s_).value()
+        << "  Max(" << s_.name() << ") = " << max(s_).value()
+        << endl;
+}
+
+
 bool Foam::functionObjects::scalarTransport::write()
 {
     return true;

src/functionObjects/solvers/scalarTransport/scalarTransport.H

@@ -30,16 +30,58 @@ Description
     - To specify the field name set the \c field entry
     - To employ the same numerical schemes as another field set
       the \c schemesField entry,
-    - A constant diffusivity may be specified with the \c D entry,
+    - The \c diffusivity entry can be set to \c none, \c constant, \c viscosity
-
+    - A constant diffusivity is specified with the \c D entry,
-    - Alternatively if a turbulence model is available a turbulent diffusivity
+    - If a momentum transport model is available and the \c viscosity
-      may be constructed from the laminar and turbulent viscosities using the
+      diffusivety option specified an effective diffusivity may be constructed
-      optional diffusivity coefficients \c alphaD and \c alphaDt (which default
+      from the laminar and turbulent viscosities using the diffusivity
-      to 1):
+      coefficients \c alphal and \c alphat:
       \verbatim
-          D = alphaD*nu + alphaDt*nut
+          D = alphal*nu + alphat*nut
       \endverbatim
 
+    Example:
+    \verbatim
+        #includeFunc scalarTransport(T, alphaD=1, alphaDt=1)
+    \endverbatim
+
+    For incompressible flow the passive scalar may optionally be solved with the
+    MULES limiter and sub-cycling or semi-implicit in order to maintain
+    boundedness, particularly if a compressive, PLIC or MPLIC convection
+    scheme is used.
+
+    Example:
+    \verbatim
+        #includeFunc scalarTransport(tracer, diffusion=none)
+
+    with scheme specification:
+        div(phi,tracer)     Gauss interfaceCompression vanLeer 1;
+
+    and solver specification:
+        tracer
+        {
+            nCorr      1;
+            nSubCycles 3;
+
+            MULESCorr       no;
+            nLimiterIter    5;
+            applyPrevCorr   yes;
+
+            solver          smoothSolver;
+            smoother        symGaussSeidel;
+            tolerance       1e-8;
+            relTol          0;
+
+            diffusion
+            {
+                solver          smoothSolver;
+                smoother        symGaussSeidel;
+                tolerance       1e-8;
+                relTol          0;
+            }
+        }
+    \endverbatim
+
 See also
     Foam::functionObjects::fvMeshFunctionObject
 
@@ -69,6 +111,19 @@ class scalarTransport
 :
     public fvMeshFunctionObject
 {
+public:
+
+        //- Enumeration defining the type of the diffusion
+        enum class diffusionType
+        {
+            none,
+            constant,
+            viscosity
+        };
+
+
+private:
+
     // Private Data
 
         //- Name of field to process
@@ -80,17 +135,20 @@ class scalarTransport
         //- Name of density field (optional)
         word rhoName_;
 
-        //- Diffusion coefficient (optional)
+        //- Diffusion type names
+        static const NamedEnum<diffusionType, 3> diffusionTypeNames_;
+
+        //- The type of diffusion
+        diffusionType diffusion_;
+
+        //- Constant diffusion coefficient (optional)
         scalar D_;
 
-        //- Flag to indicate whether a constant, uniform D_ is specified
-        bool constantD_;
-
         //- Laminar diffusion coefficient (optional)
-        scalar alphaD_;
+        scalar alphal_;
 
         //- Turbulent diffusion coefficient (optional)
-        scalar alphaDt_;
+        scalar alphat_;
 
         //- Number of corrector iterations (optional)
         int nCorr_;
@@ -101,11 +159,30 @@ class scalarTransport
         //- The scalar field
         volScalarField s_;
 
+        //- Switch for MULES limited solution
+        bool MULES_;
+
+        //- Stabilisation for normalisation of the interface normal
+        //  needed if a compressive convection scheme is used
+        const dimensionedScalar deltaN_;
+
+        //- Scalar volumetric flux
+        tmp<surfaceScalarField> tsPhi_;
+
+        //- MULES Correction
+        tmp<surfaceScalarField> tsPhiCorr0_;
+
+        //- Switch to indicate that s has been restarted for Crank-Nicolson
+        bool sRestart_;
+
 
     // Private Member Functions
 
         //- Return the diffusivity field
-        tmp<volScalarField> D(const surfaceScalarField& phi) const;
+        tmp<volScalarField> D() const;
+
+        void subCycleMULES();
+        void solveMULES();
 
 
 public:

src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H

@@ -1,13 +1,3 @@
-static const wordHashSet compressionSchemes
-{
-    "interfaceCompression",
-    "noInterfaceCompression",
-    "PLIC",
-    "PLICU",
-    "MPLIC",
-    "MPLICU"
-};
-
 static const word divAlphaName("div(phi,alpha)");
 
 const word alphaScheme(mesh.divScheme(divAlphaName)[1].wordToken());

src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2020-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,16 @@ namespace Foam
     surfaceInterpolationScheme<scalar>::
         addMeshFluxConstructorToTable<interfaceCompressionNew>
         addinterfaceCompressionScalarMeshFluxConstructorToTable_;
+
+    const wordHashSet compressionSchemes
+    {
+        "interfaceCompression",
+        "noInterfaceCompression",
+        "PLIC",
+        "PLICU",
+        "MPLIC",
+        "MPLICU"
+    };
 }
 
 

src/twoPhaseModels/twoPhaseMixture/interfaceCompression/interfaceCompression.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2020-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -143,6 +143,9 @@ public:
 };
 
 
+extern const wordHashSet compressionSchemes;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

tutorials/incompressible/pimpleFoam/RAS/flowWithOpenBoundary/system/controlDict

@@ -45,7 +45,7 @@ runTimeModifiable true;
 
 functions
 {
-    #includeFunc scalarTransport(T, alphaD=1, alphaDt=1)
+    #includeFunc scalarTransport(T, alphal=1, alphat=1)
     #includeFunc time
 }