From 16a89a6fce7fa5bed94fa96e67835a5c637b6a6b Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Fri, 25 Oct 2013 15:08:33 +0100
Subject: [PATCH] MULES: Consolidated support for LTS
---
.../multiphaseMixtureThermo.C | 1 +
.../multiphase/interFoam/LTSInterFoam/MULES.C | 83 ---
.../multiphase/interFoam/LTSInterFoam/MULES.H | 136 ----
.../interFoam/LTSInterFoam/MULESTemplates.C | 604 ------------------
.../interFoam/LTSInterFoam/Make/files | 1 -
.../interFoam/interMixingFoam/alphaEqns.H | 2 +
.../multiphaseSystem/multiphaseSystem.C | 1 +
.../multiphaseMixture/multiphaseMixture.C | 1 +
.../solvers/MULES/IMULESTemplates.C | 1 +
.../fvMatrices/solvers/MULES/MULES.C | 21 +
.../fvMatrices/solvers/MULES/MULES.H | 48 +-
.../fvMatrices/solvers/MULES/MULESTemplates.C | 126 ++--
12 files changed, 158 insertions(+), 867 deletions(-)
delete mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
delete mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
delete mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 16395505643..602affd4e00 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
+++ b/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_,
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
deleted file mode 100644
index 1720a267a0b..00000000000
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
+++ /dev/null
@@ -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
- );
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
deleted file mode 100644
index 2c2e282ae27..00000000000
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
+++ /dev/null
@@ -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
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
deleted file mode 100644
index 9185897283a..00000000000
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
+++ /dev/null
@@ -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>());
- }
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
index 205812e978d..db8b3af426e 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
@@ -1,4 +1,3 @@
LTSInterFoam.C
-MULES.C
EXE = $(FOAM_APPBIN)/LTSInterFoam
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H
index 15759adfe4a..970da9fc925 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H
+++ b/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,
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 995d07d7e45..0a276c6e6f8 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/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_,
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index 2a5d9a03b0a..20afd1f72ff 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/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_,
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C
index 14749ff9b73..542181e90dd 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C
+++ b/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,
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
index ff5994e670c..f2898b82b84 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
+++ b/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,
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
index c1d824b9503..b4253f6ee88 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
+++ b/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,
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index 53b64a5f7b5..a0d2e94daa5 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/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);
- explicitSolve(rho, psi, phiPsi, Sp, Su);
+ limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
+ 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);
- correct(rho, psi, phiCorr, Sp, Su);
+ limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+ 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.oldTime().field()/deltaT)*psi0
+ - (rho.field()*rDeltaT - Sp.field())*psiMinn
+ + (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()*psi.internalField()/deltaT
+ - (rho.field()*rDeltaT - Sp.field())*psiMinn
+ + 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,
--
GitLab