diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C index 1639550564..602affd4e0 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 1720a267a0..0000000000 --- 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 . - -\*---------------------------------------------------------------------------*/ - -#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 2c2e282ae2..0000000000 --- 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 . - -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 -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 -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 -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 9185897283..0000000000 --- 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 . - -\*---------------------------------------------------------------------------*/ - -#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 -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(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("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 -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 CDs(mesh); - upwind UDs(mesh, phi); - //fv::uncorrectedSnGrad snGrads(mesh); - - fvScalarMatrix psiConvectionDiffusion - ( - fvm::ddt(rho, psi) - + fv::gaussConvectionScheme(mesh, phi, UDs).fvmDiv(phi, psi) - //- fv::gaussLaplacianScheme(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("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 -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 tVsc = mesh.Vsc(); - const scalarField& V = tVsc(); - - const volScalarField& rDeltaT = - mesh.objectRegistry::lookupObject("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 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 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()); - } -} - - -// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files index 205812e978..db8b3af426 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 15759adfe4..970da9fc92 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 995d07d7e4..0a276c6e6f 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 2a5d9a03b0..20afd1f72f 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 14749ff9b7..542181e90d 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 ff5994e670..f2898b82b8 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 c1d824b950..b4253f6ee8 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H @@ -60,6 +60,17 @@ namespace MULES // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +template +void explicitSolve +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su +); + template void explicitSolve ( @@ -92,10 +103,33 @@ void explicitSolve const scalar psiMin ); - template +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 void correct ( + const RdeltaTType& rDeltaT, const RhoType& rho, volScalarField& psi, const surfaceScalarField& phiCorr, @@ -124,10 +158,11 @@ void correct ); -template +template 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 +template void limit ( + const RdeltaTType& rDeltaT, const RhoType& rho, const volScalarField& psi, const surfaceScalarField& phi, @@ -155,10 +191,11 @@ void limit ); -template +template 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 +template 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 53b64a5f7b..a0d2e94daa 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C @@ -32,9 +32,10 @@ License // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -template +template 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 +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 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 +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("rSubDeltaT"); + + psi.correctBoundaryConditions(); + limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false); + explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su); +} + + +template 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 +template 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 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 +template 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 +template 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 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 +template 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,