MULES: Consolidated support for LTS

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C

@@ -986,6 +986,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
 
         MULES::limit
         (
+            1.0/mesh_.time().deltaT().value(),
             geometricOneField(),
             alpha,
             phi_,

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

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

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

@@ -1,136 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Global
-    MULES
-
-Description
-    MULES: Multidimensional universal limiter with explicit solution.
-
-    Solve a convective-only transport equation using an explicit universal
-    multi-dimensional limiter.
-
-    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
-    MULES.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef MULES_H
-#define MULES_H
-
-#include "volFields.H"
-#include "surfaceFieldsFwd.H"
-#include "primitiveFieldsFwd.H"
-#include "zeroField.H"
-#include "geometricOneField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace MULES
-{
-
-template<class RhoType, class SpType, class SuType>
-void explicitLTSSolve
-(
-    const RhoType& rho,
-    volScalarField& psi,
-    const surfaceScalarField& phiBD,
-    surfaceScalarField& phiPsi,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-void explicitLTSSolve
-(
-    volScalarField& psi,
-    const surfaceScalarField& phiBD,
-    surfaceScalarField& phiPsi,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-template<class RhoType, class SpType, class SuType>
-void implicitSolve
-(
-    const RhoType& rho,
-    volScalarField& gamma,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiCorr,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-void implicitSolve
-(
-    volScalarField& gamma,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiCorr,
-    const scalar psiMax,
-    const scalar psiMin
-);
-
-template<class RhoType, class SpType, class SuType>
-void limiter
-(
-    scalarField& allLambda,
-    const RhoType& rho,
-    const volScalarField& psi,
-    const surfaceScalarField& phiBD,
-    const 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 "MULESTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

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

@@ -1,604 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "MULES.H"
-#include "upwind.H"
-#include "uncorrectedSnGrad.H"
-#include "gaussConvectionScheme.H"
-#include "gaussLaplacianScheme.H"
-#include "uncorrectedSnGrad.H"
-#include "surfaceInterpolate.H"
-#include "fvcSurfaceIntegrate.H"
-#include "slicedSurfaceFields.H"
-#include "syncTools.H"
-
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::explicitLTSSolve
-(
-    const RhoType& rho,
-    volScalarField& psi,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiPsi,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    Info<< "MULES: Solving for " << psi.name() << endl;
-
-    const fvMesh& mesh = psi.mesh();
-    psi.correctBoundaryConditions();
-
-    surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
-
-    surfaceScalarField& phiCorr = phiPsi;
-    phiCorr -= phiBD;
-
-    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
-    );
-
-    limiter
-    (
-        allLambda,
-        rho,
-        psi,
-        phiBD,
-        phiCorr,
-        Sp.field(),
-        Su.field(),
-        psiMax,
-        psiMin,
-        3
-    );
-
-    phiPsi = phiBD + lambda*phiCorr;
-
-    scalarField& psiIf = psi;
-    const scalarField& psi0 = psi.oldTime();
-
-    const volScalarField& rDeltaT =
-        mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
-
-    psiIf = 0.0;
-    fvc::surfaceIntegrate(psiIf, phiPsi);
-
-    if (mesh.moving())
-    {
-        psiIf =
-        (
-            mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc()
-          + Su.field()
-          - psiIf
-        )/(rho*rDeltaT - Sp.field());
-    }
-    else
-    {
-        psiIf =
-        (
-            rho.oldTime()*psi0*rDeltaT
-          + Su.field()
-          - psiIf
-        )/(rho*rDeltaT - Sp.field());
-    }
-
-    psi.correctBoundaryConditions();
-}
-
-
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::implicitSolve
-(
-    const RhoType& rho,
-    volScalarField& psi,
-    const surfaceScalarField& phi,
-    surfaceScalarField& phiPsi,
-    const SpType& Sp,
-    const SuType& Su,
-    const scalar psiMax,
-    const scalar psiMin
-)
-{
-    const fvMesh& mesh = psi.mesh();
-
-    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
-
-    label maxIter
-    (
-        readLabel(MULEScontrols.lookup("maxIter"))
-    );
-
-    label nLimiterIter
-    (
-        readLabel(MULEScontrols.lookup("nLimiterIter"))
-    );
-
-    scalar maxUnboundedness
-    (
-        readScalar(MULEScontrols.lookup("maxUnboundedness"))
-    );
-
-    scalar CoCoeff
-    (
-        readScalar(MULEScontrols.lookup("CoCoeff"))
-    );
-
-    scalarField allCoLambda(mesh.nFaces());
-
-    {
-        surfaceScalarField Cof
-        (
-            mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
-           *mag(phi)/mesh.magSf()
-        );
-
-        slicedSurfaceScalarField CoLambda
-        (
-            IOobject
-            (
-                "CoLambda",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimless,
-            allCoLambda,
-            false   // Use slices for the couples
-        );
-
-        CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
-    }
-
-    scalarField allLambda(allCoLambda);
-    //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
-    );
-
-    linear<scalar> CDs(mesh);
-    upwind<scalar> UDs(mesh, phi);
-    //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
-
-    fvScalarMatrix psiConvectionDiffusion
-    (
-        fvm::ddt(rho, psi)
-      + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
-        //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
-        //.fvmLaplacian(Dpsif, psi)
-      - fvm::Sp(Sp, psi)
-      - Su
-    );
-
-    surfaceScalarField phiBD(psiConvectionDiffusion.flux());
-
-    surfaceScalarField& phiCorr = phiPsi;
-    phiCorr -= phiBD;
-
-    for (label i=0; i<maxIter; i++)
-    {
-        if (i != 0 && i < 4)
-        {
-            allLambda = allCoLambda;
-        }
-
-        limiter
-        (
-            allLambda,
-            rho,
-            psi,
-            phiBD,
-            phiCorr,
-            Sp.field(),
-            Su.field(),
-            psiMax,
-            psiMin,
-            nLimiterIter
-        );
-
-        solve
-        (
-            psiConvectionDiffusion + fvc::div(lambda*phiCorr),
-            MULEScontrols
-        );
-
-        scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
-        scalar minPsi = gMin(psi.internalField());
-
-        scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
-
-        if (unboundedness < maxUnboundedness)
-        {
-            break;
-        }
-        else
-        {
-            Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
-                << " min(" << psi.name() << ") = " << minPsi << endl;
-
-            phiBD = psiConvectionDiffusion.flux();
-
-            /*
-            word gammaScheme("div(phi,gamma)");
-            word gammarScheme("div(phirb,gamma)");
-
-            const surfaceScalarField& phir =
-                mesh.lookupObject<surfaceScalarField>("phir");
-
-            phiCorr =
-                fvc::flux
-                (
-                    phi,
-                    psi,
-                    gammaScheme
-                )
-              + fvc::flux
-                (
-                    -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
-                    psi,
-                    gammarScheme
-                )
-                - phiBD;
-            */
-        }
-    }
-
-    phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
-}
-
-
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::limiter
-(
-    scalarField& allLambda,
-    const RhoType& rho,
-    const volScalarField& psi,
-    const surfaceScalarField& phiBD,
-    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 scalarField& psi0 = psi.oldTime();
-
-    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 volScalarField& rDeltaT =
-        mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
-
-    const scalarField& phiBDIf = phiBD;
-    const surfaceScalarField::GeometricBoundaryField& phiBDBf =
-        phiBD.boundaryField();
-
-    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 sumPhiBD(psiIf.size(), 0.0);
-
-    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]);
-
-        sumPhiBD[own] += phiBDIf[facei];
-        sumPhiBD[nei] -= phiBDIf[facei];
-
-        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& phiBDPf = phiBDBf[patchi];
-        const scalarField& phiCorrPf = phiCorrBf[patchi];
-
-        const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
-
-        if (psiPf.coupled())
-        {
-            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];
-
-            sumPhiBD[pfCelli] += phiBDPf[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);
-
-    if (mesh.moving())
-    {
-        tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
-
-        psiMaxn =
-            V*((rho*rDeltaT - Sp)*psiMaxn - Su)
-          - (V0()*rDeltaT)*rho.oldTime()*psi0
-          + sumPhiBD;
-
-        psiMinn =
-            V*(Su - (rho*rDeltaT - Sp)*psiMinn)
-          + (V0*rDeltaT)*rho.oldTime()*psi0
-          - sumPhiBD;
-    }
-    else
-    {
-        psiMaxn =
-            V*((rho*rDeltaT - Sp)*psiMaxn - (rho.oldTime()*rDeltaT)*psi0 - Su)
-          + sumPhiBD;
-
-        psiMinn =
-            V*((rho*rDeltaT)*psi0 - (rho.oldTime()*rDeltaT - Sp)*psiMinn + Su)
-          - sumPhiBD;
-    }
-
-    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],
-                    1.0), 0.0
-                );
-
-            mSumlPhim[celli] =
-                max(min
-                (
-                    (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
-                    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 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]);
-                }
-            }
-        }
-
-        syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
-    }
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files

@@ -1,4 +1,3 @@
 LTSInterFoam.C
-MULES.C
 
 EXE = $(FOAM_APPBIN)/LTSInterFoam

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H

@@ -75,6 +75,7 @@
         MULES::limiter
         (
             allLambda,
+            1.0/runTime.deltaT().value(),
             geometricOneField(),
             alpha1,
             phiAlpha1BD,
@@ -116,6 +117,7 @@
         MULES::limiter
         (
             allLambda,
+            1.0/runTime.deltaT().value(),
             geometricOneField(),
             alpha2,
             phiAlpha2BD,

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -148,6 +148,7 @@ void Foam::multiphaseSystem::solveAlphas()
 
         MULES::limit
         (
+            1.0/mesh_.time().deltaT().value(),
             geometricOneField(),
             phase1,
             phi_,

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -598,6 +598,7 @@ void Foam::multiphaseMixture::solveAlphas
 
         MULES::limit
         (
+            1.0/mesh_.time().deltaT().value(),
             geometricOneField(),
             alpha,
             phi_,

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

@@ -179,6 +179,7 @@ void Foam::MULES::implicitSolve
         limiter
         (
             allLambda,
+            1.0/mesh.time().deltaTValue(),
             rho,
             psi,
             phiBD,

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

@@ -48,6 +48,27 @@ void Foam::MULES::explicitSolve
 }
 
 
+void Foam::MULES::explicitLTSSolve
+(
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    explicitLTSSolve
+    (
+        geometricOneField(),
+        psi,
+        phi,
+        phiPsi,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
+
+
 void Foam::MULES::correct
 (
     volScalarField& psi,

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

@@ -60,6 +60,17 @@ namespace MULES
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void explicitSolve
+(
+    const RdeltaTType& rDeltaT,
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su
+);
+
 template<class RhoType, class SpType, class SuType>
 void explicitSolve
 (
@@ -92,10 +103,33 @@ void explicitSolve
     const scalar psiMin
 );
 
-
 template<class RhoType, class SpType, class SuType>
+void explicitLTSSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+void explicitLTSSolve
+(
+    volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void correct
 (
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
     const surfaceScalarField& phiCorr,
@@ -124,10 +158,11 @@ void correct
 );
 
 
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void limiter
 (
     scalarField& allLambda,
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
     const surfaceScalarField& phiBD,
@@ -139,9 +174,10 @@ void limiter
     const label nLimiterIter
 );
 
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void limit
 (
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
     const surfaceScalarField& phi,
@@ -155,10 +191,11 @@ void limit
 );
 
 
-template<class RhoType, class SpType, class SuType>
+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,
@@ -169,9 +206,10 @@ void limiterCorr
     const label nLimiterIter
 );
 
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void limitCorr
 (
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
     surfaceScalarField& phiCorr,

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

@@ -32,9 +32,10 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void Foam::MULES::explicitSolve
 (
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
     const surfaceScalarField& phiPsi,
@@ -48,7 +49,6 @@ void Foam::MULES::explicitSolve
 
     scalarField& psiIf = psi;
     const scalarField& psi0 = psi.oldTime();
-    const scalar deltaT = mesh.time().deltaTValue();
 
     psiIf = 0.0;
     fvc::surfaceIntegrate(psiIf, phiPsi);
@@ -58,25 +58,41 @@ void Foam::MULES::explicitSolve
         psiIf =
         (
             mesh.Vsc0()().field()*rho.oldTime().field()
-           *psi0/(deltaT*mesh.Vsc()().field())
+           *psi0*rDeltaT/mesh.Vsc()().field()
           + Su.field()
           - psiIf
-        )/(rho.field()/deltaT - Sp.field());
+        )/(rho.field()*rDeltaT - Sp.field());
     }
     else
     {
         psiIf =
         (
-            rho.oldTime().field()*psi0/deltaT
+            rho.oldTime().field()*psi0*rDeltaT
           + Su.field()
           - psiIf
-        )/(rho.field()/deltaT - Sp.field());
+        )/(rho.field()*rDeltaT - Sp.field());
     }
 
     psi.correctBoundaryConditions();
 }
 
 
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::explicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su
+)
+{
+    const fvMesh& mesh = psi.mesh();
+    const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
+    explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
+}
+
+
 template<class RhoType, class SpType, class SuType>
 void Foam::MULES::explicitSolve
 (
@@ -90,15 +106,42 @@ void Foam::MULES::explicitSolve
     const scalar psiMin
 )
 {
+    const fvMesh& mesh = psi.mesh();
+    const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
     psi.correctBoundaryConditions();
-    limit(rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
+    limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
-    explicitSolve(rho, psi, phiPsi, Sp, Su);
+    explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
 }
 
 
 template<class RhoType, class SpType, class SuType>
+void Foam::MULES::explicitLTSSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    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");
+
+    psi.correctBoundaryConditions();
+    limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
+    explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
+}
+
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void Foam::MULES::correct
 (
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     volScalarField& psi,
     const surfaceScalarField& phiCorr,
@@ -110,8 +153,6 @@ void Foam::MULES::correct
 
     const fvMesh& mesh = psi.mesh();
 
-    const scalar deltaT = mesh.time().deltaTValue();
-
     scalarField psiIf(psi.size(), 0);
     fvc::surfaceIntegrate(psiIf, phiCorr);
 
@@ -119,19 +160,19 @@ void Foam::MULES::correct
     {
         psi.internalField() =
         (
-            rho.field()*psi.internalField()/deltaT
+            rho.field()*psi.internalField()*rDeltaT
           + Su.field()
           - psiIf
-        )/(rho.field()/deltaT - Sp.field());
+        )/(rho.field()*rDeltaT - Sp.field());
     }
     else
     {
         psi.internalField() =
         (
-            rho.field()*psi.internalField()/deltaT
+            rho.field()*psi.internalField()*rDeltaT
           + Su.field()
           - psiIf
-        )/(rho.field()/deltaT - Sp.field());
+        )/(rho.field()*rDeltaT - Sp.field());
     }
 
     psi.correctBoundaryConditions();
@@ -151,6 +192,7 @@ void Foam::MULES::correct
 )
 {
     const fvMesh& mesh = psi.mesh();
+    const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
 
     const dictionary& MULEScontrols = mesh.solverDict(psi.name());
 
@@ -159,15 +201,16 @@ void Foam::MULES::correct
         readLabel(MULEScontrols.lookup("nLimiterIter"))
     );
 
-    limitCorr(rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+    limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
-    correct(rho, psi, phiCorr, Sp, Su);
+    correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
 }
 
 
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void Foam::MULES::limiter
 (
     scalarField& allLambda,
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
     const surfaceScalarField& phiBD,
@@ -190,7 +233,6 @@ void Foam::MULES::limiter
     const labelUList& neighb = mesh.neighbour();
     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
     const scalarField& V = tVsc();
-    const scalar deltaT = mesh.time().deltaTValue();
 
     const scalarField& phiBDIf = phiBD;
     const surfaceScalarField::GeometricBoundaryField& phiBDBf =
@@ -321,19 +363,19 @@ void Foam::MULES::limiter
         psiMaxn =
             V
            *(
-               (rho.field()/deltaT - Sp.field())*psiMaxn
+               (rho.field()*rDeltaT - Sp.field())*psiMaxn
              - Su.field()
             )
-          - (V0().field()/deltaT)*rho.oldTime().field()*psi0
+          - (V0().field()*rDeltaT)*rho.oldTime().field()*psi0
           + sumPhiBD;
 
         psiMinn =
             V
            *(
                Su.field()
-             - (rho.field()/deltaT - Sp.field())*psiMinn
+             - (rho.field()*rDeltaT - Sp.field())*psiMinn
             )
-          + (V0().field()/deltaT)*rho.oldTime().field()*psi0
+          + (V0().field()*rDeltaT)*rho.oldTime().field()*psi0
           - sumPhiBD;
     }
     else
@@ -341,9 +383,9 @@ void Foam::MULES::limiter
         psiMaxn =
             V
            *(
-               (rho.field()/deltaT - Sp.field())*psiMaxn
+               (rho.field()*rDeltaT - Sp.field())*psiMaxn
              - Su.field()
-             - (rho.oldTime().field()/deltaT)*psi0
+             - (rho.oldTime().field()*rDeltaT)*psi0
             )
           + sumPhiBD;
 
@@ -351,8 +393,8 @@ void Foam::MULES::limiter
             V
            *(
                Su.field()
-             - (rho.field()/deltaT - Sp.field())*psiMinn
+             - (rho.field()*rDeltaT - Sp.field())*psiMinn
-             + (rho.oldTime().field()/deltaT)*psi0
+             + (rho.oldTime().field()*rDeltaT)*psi0
             )
           - sumPhiBD;
     }
@@ -413,14 +455,16 @@ void Foam::MULES::limiter
             sumlPhip[celli] =
                 max(min
                 (
-                    (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
+                    (sumlPhip[celli] + psiMaxn[celli])
+                   /(mSumPhim[celli] - SMALL),
                     1.0), 0.0
                 );
 
             mSumlPhim[celli] =
                 max(min
                 (
-                    (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
+                    (mSumlPhim[celli] + psiMinn[celli])
+                   /(sumPhip[celli] + SMALL),
                     1.0), 0.0
                 );
         }
@@ -514,9 +558,10 @@ void Foam::MULES::limiter
 }
 
 
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void Foam::MULES::limit
 (
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
     const surfaceScalarField& phi,
@@ -558,6 +603,7 @@ void Foam::MULES::limit
     limiter
     (
         allLambda,
+        rDeltaT,
         rho,
         psi,
         phiBD,
@@ -580,10 +626,11 @@ void Foam::MULES::limit
 }
 
 
-template<class RhoType, class SpType, class SuType>
+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,
@@ -603,7 +650,6 @@ void Foam::MULES::limiterCorr
     const labelUList& neighb = mesh.neighbour();
     tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
     const scalarField& V = tVsc();
-    const scalar deltaT = mesh.time().deltaTValue();
 
     const scalarField& phiCorrIf = phiCorr;
     const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
@@ -718,17 +764,17 @@ void Foam::MULES::limiterCorr
     psiMaxn =
         V
        *(
-           (rho.field()/deltaT - Sp.field())*psiMaxn
+           (rho.field()*rDeltaT - Sp.field())*psiMaxn
          - Su.field()
-         - rho.field()*psi.internalField()/deltaT
+         - rho.field()*psi.internalField()*rDeltaT
         );
 
     psiMinn =
         V
        *(
            Su.field()
-         - (rho.field()/deltaT - Sp.field())*psiMinn
+         - (rho.field()*rDeltaT - Sp.field())*psiMinn
-         + rho.field()*psi.internalField()/deltaT
+         + rho.field()*psi.internalField()*rDeltaT
         );
 
     scalarField sumlPhip(psiIf.size());
@@ -787,14 +833,16 @@ void Foam::MULES::limiterCorr
             sumlPhip[celli] =
                 max(min
                 (
-                    (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
+                    (sumlPhip[celli] + psiMaxn[celli])
+                   /(mSumPhim[celli] - SMALL),
                     1.0), 0.0
                 );
 
             mSumlPhim[celli] =
                 max(min
                 (
-                    (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
+                    (mSumlPhim[celli] + psiMinn[celli])
+                   /(sumPhip[celli] + SMALL),
                     1.0), 0.0
                 );
         }
@@ -887,9 +935,10 @@ void Foam::MULES::limiterCorr
 }
 
 
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
 void Foam::MULES::limitCorr
 (
+    const RdeltaTType& rDeltaT,
     const RhoType& rho,
     const volScalarField& psi,
     surfaceScalarField& phiCorr,
@@ -924,6 +973,7 @@ void Foam::MULES::limitCorr
     limiterCorr
     (
         allLambda,
+        rDeltaT,
         rho,
         psi,
         phiCorr,