MULES: Added support for predictor-corrector formulation to interFoam and LTSInterFoam

applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C

@@ -38,7 +38,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
+#include "CMULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "incompressibleTwoPhaseMixture.H"
@@ -79,7 +79,10 @@ int main(int argc, char *argv[])
 
         twoPhaseProperties.correct();
 
+        #define LTSSOLVE
         #include "alphaEqnSubCycle.H"
+        #undef LTSSOLVE
+
         interface.correct();
 
         turbulence->correct();

applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H

@@ -1,39 +0,0 @@
-{
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
-
-    surfaceScalarField phic(mag(phi/mesh.magSf()));
-    phic = min(interface.cAlpha()*phic, max(phic));
-    surfaceScalarField phir(phic*interface.nHatf());
-
-    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
-    {
-        surfaceScalarField phiAlpha
-        (
-            fvc::flux
-            (
-                phi,
-                alpha1,
-                alphaScheme
-            )
-          + fvc::flux
-            (
-                -fvc::flux(-phir, alpha2, alpharScheme),
-                alpha1,
-                alpharScheme
-            )
-        );
-
-        MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
-        //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
-
-        alpha2 = 1.0 - alpha1;
-        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
-    }
-
-    Info<< "Phase-1 volume fraction = "
-        << alpha1.weightedAverage(mesh.Vsc()).value()
-        << "  Min(alpha1) = " << min(alpha1).value()
-        << "  Max(alpha1) = " << max(alpha1).value()
-        << endl;
-}

applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H

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

applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C

@@ -37,7 +37,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
+#include "CMULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "incompressibleTwoPhaseMixture.H"

applications/solvers/multiphase/interFoam/alphaEqn.H

@@ -6,9 +6,39 @@
     phic = min(interface.cAlpha()*phic, max(phic));
     surfaceScalarField phir(phic*interface.nHatf());
 
+    tmp<surfaceScalarField> tphiAlpha;
+
+    if (MULESCorr)
+    {
+        fvScalarMatrix alpha1Eqn
+        (
+            #ifdef LTSSOLVE
+            fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
+            #else
+            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+            #endif
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh,
+                phi,
+                upwind<scalar>(mesh, phi)
+            ).fvmDiv(phi, alpha1)
+        );
+
+        alpha1Eqn.solve();
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).value()
+            << endl;
+
+        tphiAlpha = alpha1Eqn.flux();
+    }
+
     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
     {
-        surfaceScalarField phiAlpha
+        tmp<surfaceScalarField> tphiAlpha0
         (
             fvc::flux
             (
@@ -24,12 +54,32 @@
             )
         );
 
-        MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
+        if (MULESCorr)
+        {
+            tphiAlpha0() -= tphiAlpha();
+            #ifdef LTSSOLVE
+            MULES::LTScorrect(alpha1, tphiAlpha0(), 1, 0);
+            #else
+            MULES::correct(alpha1, tphiAlpha0(), 1, 0);
+            #endif
+            tphiAlpha() += tphiAlpha0();
+        }
+        else
+        {
+            tphiAlpha = tphiAlpha0;
+
+            #ifdef LTSSOLVE
+            MULES::explicitLTSSolve(alpha1, phi, tphiAlpha(), 1, 0);
+            #else
+            MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0);
+            #endif
+        }
 
         alpha2 = 1.0 - alpha1;
-        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
     }
 
+    rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
+
     Info<< "Phase-1 volume fraction = "
         << alpha1.weightedAverage(mesh.Vsc()).value()
         << "  Min(alpha1) = " << min(alpha1).value()

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -34,7 +34,7 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "MULES.H"
+#include "CMULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "incompressibleTwoPhaseMixture.H"

applications/solvers/multiphase/interFoam/interFoam.C

@@ -38,7 +38,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
+#include "CMULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "incompressibleTwoPhaseMixture.H"

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

@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
+#include "CMULES.H"
 #include "subCycle.H"
 #include "threePhaseMixture.H"
 #include "threePhaseInterfaceProperties.H"

applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C

@@ -39,7 +39,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
+#include "CMULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "incompressibleTwoPhaseMixture.H"

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -16,8 +16,14 @@
     {
         fvScalarMatrix alpha1Eqn
         (
-            fvm::ddt(alpha1)
+            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
-          + fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1)
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh,
+                phi,
+                upwind<scalar>(mesh, phi)
+            ).fvmDiv(phi, alpha1)
+          - fvm::Sp(divU, alpha1)
          ==
             fvm::Sp(vDotvmcAlphal, alpha1)
           + vDotcAlphal

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C

@@ -43,7 +43,7 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "MULES.H"
+#include "CMULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "phaseChangeTwoPhaseMixture.H"

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -41,7 +41,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
+#include "CMULES.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "phaseChangeTwoPhaseMixture.H"

src/finiteVolume/Make/files

@@ -220,6 +220,7 @@ fields/surfaceFields/surfaceFields.C
 fvMatrices/fvMatrices.C
 fvMatrices/fvScalarMatrix/fvScalarMatrix.C
 fvMatrices/solvers/MULES/MULES.C
+fvMatrices/solvers/MULES/CMULES.C
 fvMatrices/solvers/MULES/IMULES.C
 fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
 

src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C

@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "CMULES.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::MULES::correct
+(
+    volScalarField& psi,
+    surfaceScalarField& phiPsiCorr,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    correct
+    (
+        geometricOneField(),
+        psi,
+        phiPsiCorr,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
+
+
+void Foam::MULES::LTScorrect
+(
+    volScalarField& psi,
+    surfaceScalarField& phiPsiCorr,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    LTScorrect
+    (
+        geometricOneField(),
+        psi,
+        phiPsiCorr,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    MULES
+
+Description
+    CMULES: Multidimensional universal limiter for
+    explicit corrected implicit solution.
+
+    Solve a convective-only transport equation using an explicit universal
+    multi-dimensional limiter to correct an implicit conservative bounded
+    obtained using rigorously bounded schemes such as Euler-implicit in time
+    upwind in space.
+
+    Parameters are the variable to solve, the normal convective flux and the
+    actual explicit flux of the variable which is also used to return limited
+    flux used in the bounded-solution.
+
+SourceFiles
+    CMULES.C
+    CMULESTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CMULES_H
+#define CMULES_H
+
+#include "MULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "gaussConvectionScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace MULES
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void correct
+(
+    const RdeltaTType& rDeltaT,
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su
+);
+
+template<class RhoType, class SpType, class SuType>
+void correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+void correct
+(
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+template<class RhoType, class SpType, class SuType>
+void LTScorrect
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+void LTScorrect
+(
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void limiterCorr
+(
+    scalarField& allLambda,
+    const RdeltaTType& rDeltaT,
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+);
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void limitCorr
+(
+    const RdeltaTType& rDeltaT,
+    const RhoType& rho,
+    const volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace MULES
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "CMULESTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C

@@ -0,0 +1,493 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "CMULES.H"
+#include "fvcSurfaceIntegrate.H"
+#include "slicedSurfaceFields.H"
+#include "wedgeFvPatch.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void Foam::MULES::correct
+(
+    const RdeltaTType& rDeltaT,
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su
+)
+{
+    Info<< "MULES: Correcting " << psi.name() << endl;
+
+    const fvMesh& mesh = psi.mesh();
+
+    scalarField psiIf(psi.size(), 0);
+    fvc::surfaceIntegrate(psiIf, phiCorr);
+
+    if (mesh.moving())
+    {
+        psi.internalField() =
+        (
+            rho.field()*psi.internalField()*rDeltaT
+          + Su.field()
+          - psiIf
+        )/(rho.field()*rDeltaT - Sp.field());
+    }
+    else
+    {
+        psi.internalField() =
+        (
+            rho.field()*psi.internalField()*rDeltaT
+          + Su.field()
+          - psiIf
+        )/(rho.field()*rDeltaT - Sp.field());
+    }
+
+    psi.correctBoundaryConditions();
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::correct
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    const fvMesh& mesh = psi.mesh();
+    const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
+
+    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+    label nLimiterIter
+    (
+        readLabel(MULEScontrols.lookup("nLimiterIter"))
+    );
+
+    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::LTScorrect
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    const volScalarField& rDeltaT =
+        mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
+
+    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+    label nLimiterIter
+    (
+        readLabel(MULEScontrols.lookup("nLimiterIter"))
+    );
+
+    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+}
+
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void Foam::MULES::limiterCorr
+(
+    scalarField& allLambda,
+    const RdeltaTType& rDeltaT,
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+)
+{
+    const scalarField& psiIf = psi;
+    const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
+
+    const fvMesh& mesh = psi.mesh();
+
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighb = mesh.neighbour();
+    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
+    const scalarField& V = tVsc();
+
+    const scalarField& phiCorrIf = phiCorr;
+    const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
+        phiCorr.boundaryField();
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    scalarField& lambdaIf = lambda;
+    surfaceScalarField::GeometricBoundaryField& lambdaBf =
+        lambda.boundaryField();
+
+    scalarField psiMaxn(psiIf.size(), psiMin);
+    scalarField psiMinn(psiIf.size(), psiMax);
+
+    scalarField sumPhip(psiIf.size(), VSMALL);
+    scalarField mSumPhim(psiIf.size(), VSMALL);
+
+    forAll(phiCorrIf, facei)
+    {
+        label own = owner[facei];
+        label nei = neighb[facei];
+
+        psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
+        psiMinn[own] = min(psiMinn[own], psiIf[nei]);
+
+        psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
+        psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
+
+        scalar phiCorrf = phiCorrIf[facei];
+
+        if (phiCorrf > 0.0)
+        {
+            sumPhip[own] += phiCorrf;
+            mSumPhim[nei] += phiCorrf;
+        }
+        else
+        {
+            mSumPhim[own] -= phiCorrf;
+            sumPhip[nei] -= phiCorrf;
+        }
+    }
+
+    forAll(phiCorrBf, patchi)
+    {
+        const fvPatchScalarField& psiPf = psiBf[patchi];
+        const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+        const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+        if (psiPf.coupled())
+        {
+            const scalarField psiPNf(psiPf.patchNeighbourField());
+
+            forAll(phiCorrPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
+            }
+        }
+        else
+        {
+            forAll(phiCorrPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
+            }
+        }
+
+        forAll(phiCorrPf, pFacei)
+        {
+            label pfCelli = pFaceCells[pFacei];
+
+            scalar phiCorrf = phiCorrPf[pFacei];
+
+            if (phiCorrf > 0.0)
+            {
+                sumPhip[pfCelli] += phiCorrf;
+            }
+            else
+            {
+                mSumPhim[pfCelli] -= phiCorrf;
+            }
+        }
+    }
+
+    psiMaxn = min(psiMaxn, psiMax);
+    psiMinn = max(psiMinn, psiMin);
+
+    // scalar smooth = 0.5;
+    // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
+    // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
+
+    psiMaxn =
+        V
+       *(
+           (rho.field()*rDeltaT - Sp.field())*psiMaxn
+         - Su.field()
+         - rho.field()*psi.internalField()*rDeltaT
+        );
+
+    psiMinn =
+        V
+       *(
+           Su.field()
+         - (rho.field()*rDeltaT - Sp.field())*psiMinn
+         + rho.field()*psi.internalField()*rDeltaT
+        );
+
+    scalarField sumlPhip(psiIf.size());
+    scalarField mSumlPhim(psiIf.size());
+
+    for (int j=0; j<nLimiterIter; j++)
+    {
+        sumlPhip = 0.0;
+        mSumlPhim = 0.0;
+
+        forAll(lambdaIf, facei)
+        {
+            label own = owner[facei];
+            label nei = neighb[facei];
+
+            scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
+
+            if (lambdaPhiCorrf > 0.0)
+            {
+                sumlPhip[own] += lambdaPhiCorrf;
+                mSumlPhim[nei] += lambdaPhiCorrf;
+            }
+            else
+            {
+                mSumlPhim[own] -= lambdaPhiCorrf;
+                sumlPhip[nei] -= lambdaPhiCorrf;
+            }
+        }
+
+        forAll(lambdaBf, patchi)
+        {
+            scalarField& lambdaPf = lambdaBf[patchi];
+            const scalarField& phiCorrfPf = phiCorrBf[patchi];
+
+            const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+            forAll(lambdaPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
+
+                if (lambdaPhiCorrf > 0.0)
+                {
+                    sumlPhip[pfCelli] += lambdaPhiCorrf;
+                }
+                else
+                {
+                    mSumlPhim[pfCelli] -= lambdaPhiCorrf;
+                }
+            }
+        }
+
+        forAll(sumlPhip, celli)
+        {
+            sumlPhip[celli] =
+                max(min
+                (
+                    (sumlPhip[celli] + psiMaxn[celli])
+                   /(mSumPhim[celli] - SMALL),
+                    1.0), 0.0
+                );
+
+            mSumlPhim[celli] =
+                max(min
+                (
+                    (mSumlPhim[celli] + psiMinn[celli])
+                   /(sumPhip[celli] + SMALL),
+                    1.0), 0.0
+                );
+        }
+
+        const scalarField& lambdam = sumlPhip;
+        const scalarField& lambdap = mSumlPhim;
+
+        forAll(lambdaIf, facei)
+        {
+            if (phiCorrIf[facei] > 0.0)
+            {
+                lambdaIf[facei] = min
+                (
+                    lambdaIf[facei],
+                    min(lambdap[owner[facei]], lambdam[neighb[facei]])
+                );
+            }
+            else
+            {
+                lambdaIf[facei] = min
+                (
+                    lambdaIf[facei],
+                    min(lambdam[owner[facei]], lambdap[neighb[facei]])
+                );
+            }
+        }
+
+
+        forAll(lambdaBf, patchi)
+        {
+            fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
+            const scalarField& phiCorrfPf = phiCorrBf[patchi];
+            const fvPatchScalarField& psiPf = psiBf[patchi];
+
+            if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
+            {
+                lambdaPf = 0;
+            }
+            else if (psiPf.coupled())
+            {
+                const labelList& pFaceCells =
+                    mesh.boundary()[patchi].faceCells();
+
+                forAll(lambdaPf, pFacei)
+                {
+                    label pfCelli = pFaceCells[pFacei];
+
+                    if (phiCorrfPf[pFacei] > 0.0)
+                    {
+                        lambdaPf[pFacei] =
+                            min(lambdaPf[pFacei], lambdap[pfCelli]);
+                    }
+                    else
+                    {
+                        lambdaPf[pFacei] =
+                            min(lambdaPf[pFacei], lambdam[pfCelli]);
+                    }
+                }
+            }
+            else
+            {
+                const labelList& pFaceCells =
+                    mesh.boundary()[patchi].faceCells();
+                const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+                forAll(lambdaPf, pFacei)
+                {
+                    // Limit outlet faces only
+                    if (phiCorrPf[pFacei] > SMALL*SMALL)
+                    {
+                        label pfCelli = pFaceCells[pFacei];
+
+                        if (phiCorrfPf[pFacei] > 0.0)
+                        {
+                            lambdaPf[pFacei] =
+                                min(lambdaPf[pFacei], lambdap[pfCelli]);
+                        }
+                        else
+                        {
+                            lambdaPf[pFacei] =
+                                min(lambdaPf[pFacei], lambdam[pfCelli]);
+                        }
+                    }
+                }
+            }
+        }
+
+        syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
+    }
+}
+
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void Foam::MULES::limitCorr
+(
+    const RdeltaTType& rDeltaT,
+    const RhoType& rho,
+    const volScalarField& psi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    scalarField allLambda(mesh.nFaces(), 1.0);
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    limiterCorr
+    (
+        allLambda,
+        rDeltaT,
+        rho,
+        psi,
+        phiCorr,
+        Sp,
+        Su,
+        psiMax,
+        psiMin,
+        nLimiterIter
+    );
+
+    phiCorr *= lambda;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H

@@ -25,7 +25,14 @@ Global
     IMULES
 
 Description
-    IMULES: Multidimensional universal limiter with implicit solution.
+    IMULES: Multidimensional universal limiter for implicit solution.
+
+    Solve a convective-only transport equation using an explicit universal
+    multi-dimensional limiter applied to an implicit formulation requiring
+    iteration to guarantee boundedness.  The number of iterations required
+    to obtain boundedness increases with the Courant number of the simulation.
+
+    It may be more efficient to use CMULES.
 
 SourceFiles
     IMULES.C

src/finiteVolume/fvMatrices/solvers/MULES/MULES.C

@@ -69,25 +69,6 @@ void Foam::MULES::explicitLTSSolve
 }
 
 
-void Foam::MULES::correct
-(
-    volScalarField& psi,
-    surfaceScalarField& phiPsiCorr,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    correct
-    (
-        geometricOneField(),
-        psi,
-        phiPsiCorr,
-        zeroField(), zeroField(),
-        psiMax, psiMin
-    );
-}
-
-
 void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
 {
     forAll(phiPsiCorrs[0], facei)

src/finiteVolume/fvMatrices/solvers/MULES/MULES.H

@@ -25,7 +25,7 @@ Global
     MULES
 
 Description
-    MULES: Multidimensional universal limiter with explicit solution.
+    MULES: Multidimensional universal limiter for explicit solution.
 
     Solve a convective-only transport equation using an explicit universal
     multi-dimensional limiter.
@@ -126,38 +126,6 @@ void explicitLTSSolve
 );
 
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
-void correct
-(
-    const RdeltaTType& rDeltaT,
-    const RhoType& rho,
-    volScalarField& psi,
-    const surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su
-);
-
-template<class RhoType, class SpType, class SuType>
-void correct
-(
-    const RhoType& rho,
-    volScalarField& psi,
-    surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-void correct
-(
-    volScalarField& psi,
-    surfaceScalarField& phiCorr,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-
 template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void limiter
 (
@@ -191,41 +159,12 @@ void limit
 );
 
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
-void limiterCorr
-(
-    scalarField& allLambda,
-    const RdeltaTType& rDeltaT,
-    const RhoType& rho,
-    const volScalarField& psi,
-    const surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin,
-    const label nLimiterIter
-);
-
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
-void limitCorr
-(
-    const RdeltaTType& rDeltaT,
-    const RhoType& rho,
-    const volScalarField& psi,
-    surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin,
-    const label nLimiterIter
-);
-
-
 void limitSum(UPtrList<scalarField>& phiPsiCorrs);
 
 template<class SurfaceScalarFieldList>
 void limitSum(SurfaceScalarFieldList& phiPsiCorrs);
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace MULES

src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C

@@ -138,74 +138,6 @@ void Foam::MULES::explicitLTSSolve
 }
 
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
-void Foam::MULES::correct
-(
-    const RdeltaTType& rDeltaT,
-    const RhoType& rho,
-    volScalarField& psi,
-    const surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su
-)
-{
-    Info<< "MULES: Correcting " << psi.name() << endl;
-
-    const fvMesh& mesh = psi.mesh();
-
-    scalarField psiIf(psi.size(), 0);
-    fvc::surfaceIntegrate(psiIf, phiCorr);
-
-    if (mesh.moving())
-    {
-        psi.internalField() =
-        (
-            rho.field()*psi.internalField()*rDeltaT
-          + Su.field()
-          - psiIf
-        )/(rho.field()*rDeltaT - Sp.field());
-    }
-    else
-    {
-        psi.internalField() =
-        (
-            rho.field()*psi.internalField()*rDeltaT
-          + Su.field()
-          - psiIf
-        )/(rho.field()*rDeltaT - Sp.field());
-    }
-
-    psi.correctBoundaryConditions();
-}
-
-
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::correct
-(
-    const RhoType& rho,
-    volScalarField& psi,
-    surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    const fvMesh& mesh = psi.mesh();
-    const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
-
-    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
-
-    label nLimiterIter
-    (
-        readLabel(MULEScontrols.lookup("nLimiterIter"))
-    );
-
-    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
-    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
-}
-
-
 template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void Foam::MULES::limiter
 (
@@ -626,368 +558,6 @@ void Foam::MULES::limit
 }
 
 
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
-void Foam::MULES::limiterCorr
-(
-    scalarField& allLambda,
-    const RdeltaTType& rDeltaT,
-    const RhoType& rho,
-    const volScalarField& psi,
-    const surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin,
-    const label nLimiterIter
-)
-{
-    const scalarField& psiIf = psi;
-    const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
-
-    const fvMesh& mesh = psi.mesh();
-
-    const labelUList& owner = mesh.owner();
-    const labelUList& neighb = mesh.neighbour();
-    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
-    const scalarField& V = tVsc();
-
-    const scalarField& phiCorrIf = phiCorr;
-    const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
-        phiCorr.boundaryField();
-
-    slicedSurfaceScalarField lambda
-    (
-        IOobject
-        (
-            "lambda",
-            mesh.time().timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh,
-        dimless,
-        allLambda,
-        false   // Use slices for the couples
-    );
-
-    scalarField& lambdaIf = lambda;
-    surfaceScalarField::GeometricBoundaryField& lambdaBf =
-        lambda.boundaryField();
-
-    scalarField psiMaxn(psiIf.size(), psiMin);
-    scalarField psiMinn(psiIf.size(), psiMax);
-
-    scalarField sumPhip(psiIf.size(), VSMALL);
-    scalarField mSumPhim(psiIf.size(), VSMALL);
-
-    forAll(phiCorrIf, facei)
-    {
-        label own = owner[facei];
-        label nei = neighb[facei];
-
-        psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
-        psiMinn[own] = min(psiMinn[own], psiIf[nei]);
-
-        psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
-        psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
-
-        scalar phiCorrf = phiCorrIf[facei];
-
-        if (phiCorrf > 0.0)
-        {
-            sumPhip[own] += phiCorrf;
-            mSumPhim[nei] += phiCorrf;
-        }
-        else
-        {
-            mSumPhim[own] -= phiCorrf;
-            sumPhip[nei] -= phiCorrf;
-        }
-    }
-
-    forAll(phiCorrBf, patchi)
-    {
-        const fvPatchScalarField& psiPf = psiBf[patchi];
-        const scalarField& phiCorrPf = phiCorrBf[patchi];
-
-        const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
-
-        if (psiPf.coupled())
-        {
-            const scalarField psiPNf(psiPf.patchNeighbourField());
-
-            forAll(phiCorrPf, pFacei)
-            {
-                label pfCelli = pFaceCells[pFacei];
-
-                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
-                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
-            }
-        }
-        else
-        {
-            forAll(phiCorrPf, pFacei)
-            {
-                label pfCelli = pFaceCells[pFacei];
-
-                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
-                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
-            }
-        }
-
-        forAll(phiCorrPf, pFacei)
-        {
-            label pfCelli = pFaceCells[pFacei];
-
-            scalar phiCorrf = phiCorrPf[pFacei];
-
-            if (phiCorrf > 0.0)
-            {
-                sumPhip[pfCelli] += phiCorrf;
-            }
-            else
-            {
-                mSumPhim[pfCelli] -= phiCorrf;
-            }
-        }
-    }
-
-    psiMaxn = min(psiMaxn, psiMax);
-    psiMinn = max(psiMinn, psiMin);
-
-    // scalar smooth = 0.5;
-    // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
-    // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
-
-    psiMaxn =
-        V
-       *(
-           (rho.field()*rDeltaT - Sp.field())*psiMaxn
-         - Su.field()
-         - rho.field()*psi.internalField()*rDeltaT
-        );
-
-    psiMinn =
-        V
-       *(
-           Su.field()
-         - (rho.field()*rDeltaT - Sp.field())*psiMinn
-         + rho.field()*psi.internalField()*rDeltaT
-        );
-
-    scalarField sumlPhip(psiIf.size());
-    scalarField mSumlPhim(psiIf.size());
-
-    for (int j=0; j<nLimiterIter; j++)
-    {
-        sumlPhip = 0.0;
-        mSumlPhim = 0.0;
-
-        forAll(lambdaIf, facei)
-        {
-            label own = owner[facei];
-            label nei = neighb[facei];
-
-            scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
-
-            if (lambdaPhiCorrf > 0.0)
-            {
-                sumlPhip[own] += lambdaPhiCorrf;
-                mSumlPhim[nei] += lambdaPhiCorrf;
-            }
-            else
-            {
-                mSumlPhim[own] -= lambdaPhiCorrf;
-                sumlPhip[nei] -= lambdaPhiCorrf;
-            }
-        }
-
-        forAll(lambdaBf, patchi)
-        {
-            scalarField& lambdaPf = lambdaBf[patchi];
-            const scalarField& phiCorrfPf = phiCorrBf[patchi];
-
-            const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
-
-            forAll(lambdaPf, pFacei)
-            {
-                label pfCelli = pFaceCells[pFacei];
-
-                scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
-
-                if (lambdaPhiCorrf > 0.0)
-                {
-                    sumlPhip[pfCelli] += lambdaPhiCorrf;
-                }
-                else
-                {
-                    mSumlPhim[pfCelli] -= lambdaPhiCorrf;
-                }
-            }
-        }
-
-        forAll(sumlPhip, celli)
-        {
-            sumlPhip[celli] =
-                max(min
-                (
-                    (sumlPhip[celli] + psiMaxn[celli])
-                   /(mSumPhim[celli] - SMALL),
-                    1.0), 0.0
-                );
-
-            mSumlPhim[celli] =
-                max(min
-                (
-                    (mSumlPhim[celli] + psiMinn[celli])
-                   /(sumPhip[celli] + SMALL),
-                    1.0), 0.0
-                );
-        }
-
-        const scalarField& lambdam = sumlPhip;
-        const scalarField& lambdap = mSumlPhim;
-
-        forAll(lambdaIf, facei)
-        {
-            if (phiCorrIf[facei] > 0.0)
-            {
-                lambdaIf[facei] = min
-                (
-                    lambdaIf[facei],
-                    min(lambdap[owner[facei]], lambdam[neighb[facei]])
-                );
-            }
-            else
-            {
-                lambdaIf[facei] = min
-                (
-                    lambdaIf[facei],
-                    min(lambdam[owner[facei]], lambdap[neighb[facei]])
-                );
-            }
-        }
-
-
-        forAll(lambdaBf, patchi)
-        {
-            fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
-            const scalarField& phiCorrfPf = phiCorrBf[patchi];
-            const fvPatchScalarField& psiPf = psiBf[patchi];
-
-            if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
-            {
-                lambdaPf = 0;
-            }
-            else if (psiPf.coupled())
-            {
-                const labelList& pFaceCells =
-                    mesh.boundary()[patchi].faceCells();
-
-                forAll(lambdaPf, pFacei)
-                {
-                    label pfCelli = pFaceCells[pFacei];
-
-                    if (phiCorrfPf[pFacei] > 0.0)
-                    {
-                        lambdaPf[pFacei] =
-                            min(lambdaPf[pFacei], lambdap[pfCelli]);
-                    }
-                    else
-                    {
-                        lambdaPf[pFacei] =
-                            min(lambdaPf[pFacei], lambdam[pfCelli]);
-                    }
-                }
-            }
-            else
-            {
-                const labelList& pFaceCells =
-                    mesh.boundary()[patchi].faceCells();
-                const scalarField& phiCorrPf = phiCorrBf[patchi];
-
-                forAll(lambdaPf, pFacei)
-                {
-                    // Limit outlet faces only
-                    if (phiCorrPf[pFacei] > SMALL*SMALL)
-                    {
-                        label pfCelli = pFaceCells[pFacei];
-
-                        if (phiCorrfPf[pFacei] > 0.0)
-                        {
-                            lambdaPf[pFacei] =
-                                min(lambdaPf[pFacei], lambdap[pfCelli]);
-                        }
-                        else
-                        {
-                            lambdaPf[pFacei] =
-                                min(lambdaPf[pFacei], lambdam[pfCelli]);
-                        }
-                    }
-                }
-            }
-        }
-
-        syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
-    }
-}
-
-
-template<class RdeltaTType, class RhoType, class SpType, class SuType>
-void Foam::MULES::limitCorr
-(
-    const RdeltaTType& rDeltaT,
-    const RhoType& rho,
-    const volScalarField& psi,
-    surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin,
-    const label nLimiterIter
-)
-{
-    const fvMesh& mesh = psi.mesh();
-
-    scalarField allLambda(mesh.nFaces(), 1.0);
-
-    slicedSurfaceScalarField lambda
-    (
-        IOobject
-        (
-            "lambda",
-            mesh.time().timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh,
-        dimless,
-        allLambda,
-        false   // Use slices for the couples
-    );
-
-    limiterCorr
-    (
-        allLambda,
-        rDeltaT,
-        rho,
-        psi,
-        phiCorr,
-        Sp,
-        Su,
-        psiMax,
-        psiMin,
-        nLimiterIter
-    );
-
-    phiCorr *= lambda;
-}
-
-
 template<class SurfaceScalarFieldList>
 void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
 {

tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict

@@ -17,7 +17,7 @@ FoamFile
 
 application     interFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
 startTime       0;
 
@@ -25,7 +25,7 @@ stopAt          endTime;
 
 endTime         100;
 
-deltaT          0.01;
+deltaT          0.1;
 
 writeControl    adjustableRunTime;
 
@@ -46,8 +46,9 @@ timePrecision   6;
 runTimeModifiable yes;
 
 adjustTimeStep  yes;
-maxCo           0.5;
+
-maxAlphaCo      0.5;
+maxCo           6;
+maxAlphaCo      6;
 maxDeltaT       1;
 
 functions

tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes

@@ -27,12 +27,13 @@ gradSchemes
 
 divSchemes
 {
-    div(rho*phi,U)  Gauss linearUpwind grad(U);
+    default             none;
-    div(phi,alpha)  Gauss vanLeer;
-    div(phirb,alpha) Gauss interfaceCompression;
 
-    div(phi,k)      Gauss upwind;
+    div(rho*phi,U)      Gauss linearUpwind grad(U);
-    div(phi,omega)  $div(phi,k);
+    div(phi,alpha)      Gauss vanLeer;
+    div(phirb,alpha)    Gauss interfaceCompression;
+
+    "div\(phi,(k|omega)\)"      Gauss upwind;
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 

tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution

@@ -20,8 +20,16 @@ solvers
     alpha.water
     {
         nAlphaCorr      1;
-        nAlphaSubCycles 3;
+        nAlphaSubCycles 1;
         cAlpha          1;
+
+        MULESCorr       yes;
+        nLimiterIter    3;
+
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
     }
 
     pcorr
@@ -30,9 +38,9 @@ solvers
         preconditioner
         {
             preconditioner  GAMG;
-            tolerance       1e-05;
+            tolerance       1e-5;
             relTol          0;
-            smoother        DICGaussSeidel;
+            smoother        GaussSeidel;
             nPreSweeps      0;
             nPostSweeps     2;
             nFinestSweeps   2;
@@ -41,31 +49,43 @@ solvers
             agglomerator    faceAreaPair;
             mergeLevels     1;
         }
-        tolerance       1e-10;
+        tolerance       1e-5;
         relTol          0;
-        maxIter         100;
+        maxIter         50;
     }
 
     p_rgh
     {
-        $pcorr;
+        solver           GAMG;
-        tolerance       1e-6;
+        tolerance        5e-9;
-        relTol          0.01;
+        relTol           0.01;
+
+        smoother         GaussSeidel;
+        nPreSweeps       0;
+        nPostSweeps      2;
+
+        cacheAgglomeration true;
+
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+
+        maxIter          50;
     };
 
     p_rghFinal
     {
         $p_rgh;
-        tolerance       1e-6;
+        tolerance       5e-9;
         relTol          0;
     }
 
     "(U|k|omega).*"
     {
         solver          smoothSolver;
-        smoother        GaussSeidel;
+        smoother        symGaussSeidel;
         nSweeps         1;
-        tolerance       1e-7;
+        tolerance       1e-6;
         relTol          0.1;
     };
 }
@@ -84,6 +104,7 @@ relaxationFactors
     }
     equations
     {
+        ".*" 1;
     }
 }
 

tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes

@@ -32,7 +32,6 @@ divSchemes
     div(phi,k)           Gauss linearUpwind grad(k);
     div(phi,alpha)       Gauss vanLeer;
     div(phirb,alpha)     Gauss interfaceCompression;
-    UD                   Gauss upwind;
     div((muEff*dev(T(grad(U))))) Gauss linear;
 }