From 16a89a6fce7fa5bed94fa96e67835a5c637b6a6b Mon Sep 17 00:00:00 2001 From: Henry Date: Fri, 25 Oct 2013 15:08:33 +0100 Subject: [PATCH 1/4] 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 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, From 45a00a1377805e4b9c7501d066317aa840f8d009 Mon Sep 17 00:00:00 2001 From: Henry Date: Fri, 25 Oct 2013 15:08:54 +0100 Subject: [PATCH 2/4] interstitialInletVelocity: corrected mapping for decomposition and reconstruction --- ...erstitialInletVelocityFvPatchVectorField.C | 39 +++++++++++++++---- ...erstitialInletVelocityFvPatchVectorField.H | 22 ++++++++++- 2 files changed, 52 insertions(+), 9 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C index 745918403b..e5f7a049f4 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C @@ -38,7 +38,7 @@ interstitialInletVelocityFvPatchVectorField const DimensionedField& iF ) : - fixedValueFvPatchField(p, iF), + fixedValueFvPatchVectorField(p, iF), inletVelocity_(p.size(), vector::zero), alphaName_("alpha") {} @@ -53,8 +53,8 @@ interstitialInletVelocityFvPatchVectorField const fvPatchFieldMapper& mapper ) : - fixedValueFvPatchField(ptf, p, iF, mapper), - inletVelocity_(ptf.inletVelocity_), + fixedValueFvPatchVectorField(ptf, p, iF, mapper), + inletVelocity_(ptf.inletVelocity_, mapper), alphaName_(ptf.alphaName_) {} @@ -67,7 +67,7 @@ interstitialInletVelocityFvPatchVectorField const dictionary& dict ) : - fixedValueFvPatchField(p, iF, dict), + fixedValueFvPatchVectorField(p, iF, dict), inletVelocity_("inletVelocity", dict, p.size()), alphaName_(dict.lookupOrDefault("alpha", "alpha")) {} @@ -79,7 +79,7 @@ interstitialInletVelocityFvPatchVectorField const interstitialInletVelocityFvPatchVectorField& ptf ) : - fixedValueFvPatchField(ptf), + fixedValueFvPatchVectorField(ptf), inletVelocity_(ptf.inletVelocity_), alphaName_(ptf.alphaName_) {} @@ -92,7 +92,7 @@ interstitialInletVelocityFvPatchVectorField const DimensionedField& iF ) : - fixedValueFvPatchField(ptf, iF), + fixedValueFvPatchVectorField(ptf, iF), inletVelocity_(ptf.inletVelocity_), alphaName_(ptf.alphaName_) {} @@ -100,6 +100,31 @@ interstitialInletVelocityFvPatchVectorField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +void Foam::interstitialInletVelocityFvPatchVectorField::autoMap +( + const fvPatchFieldMapper& m +) +{ + fixedValueFvPatchVectorField::autoMap(m); + inletVelocity_.autoMap(m); +} + + +void Foam::interstitialInletVelocityFvPatchVectorField::rmap +( + const fvPatchVectorField& ptf, + const labelList& addr +) +{ + fixedValueFvPatchVectorField::rmap(ptf, addr); + + const interstitialInletVelocityFvPatchVectorField& tiptf = + refCast(ptf); + + inletVelocity_.rmap(tiptf.inletVelocity_, addr); +} + + void Foam::interstitialInletVelocityFvPatchVectorField::updateCoeffs() { if (updated()) @@ -111,7 +136,7 @@ void Foam::interstitialInletVelocityFvPatchVectorField::updateCoeffs() patch().lookupPatchField(alphaName_); operator==(inletVelocity_/alphap); - fixedValueFvPatchField::updateCoeffs(); + fixedValueFvPatchVectorField::updateCoeffs(); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.H index 1952b82633..6affdefada 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.H @@ -141,8 +141,26 @@ public: // Member functions - //- Update the coefficients associated with the patch field - virtual void updateCoeffs(); + // Mapping functions + + //- Map (and resize as needed) from self given a mapping object + virtual void autoMap + ( + const fvPatchFieldMapper& + ); + + //- Reverse map the given fvPatchField onto this fvPatchField + virtual void rmap + ( + const fvPatchVectorField&, + const labelList& + ); + + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); //- Write virtual void write(Ostream&) const; From 7c1bc90fc5433eff7b3758063e0d798227b96fdc Mon Sep 17 00:00:00 2001 From: Henry Date: Sun, 27 Oct 2013 20:59:31 +0000 Subject: [PATCH 3/4] ddtSchemes: Added missing specialisation --- .../CoEulerDdtScheme/CoEulerDdtScheme.H | 8 +++++++- .../CrankNicolsonDdtScheme.H | 8 +++++++- .../ddtSchemes/EulerDdtScheme/EulerDdtScheme.H | 8 +++++++- .../ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H | 8 +++++++- .../backwardDdtScheme/backwardDdtScheme.H | 8 +++++++- .../boundedDdtScheme/boundedDdtScheme.H | 8 +++++++- .../localEulerDdtScheme/localEulerDdtScheme.H | 15 ++++++++++++++- .../steadyStateDdtScheme/steadyStateDdtScheme.H | 8 +++++++- 8 files changed, 63 insertions(+), 8 deletions(-) diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H index 58a6d067a9..204a793587 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H @@ -198,7 +198,6 @@ tmp CoEulerDdtScheme::fvcDdtUfCorr const GeometricField& Uf ); - template<> tmp CoEulerDdtScheme::fvcDdtPhiCorr ( @@ -206,6 +205,13 @@ tmp CoEulerDdtScheme::fvcDdtPhiCorr const surfaceScalarField& phi ); +template<> +tmp CoEulerDdtScheme::fvcDdtUfCorr +( + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& Uf +); template<> tmp CoEulerDdtScheme::fvcDdtPhiCorr diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H index 41fbbd0c28..94c34b8ea8 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H @@ -270,7 +270,6 @@ tmp CrankNicolsonDdtScheme::fvcDdtUfCorr const GeometricField& Uf ); - template<> tmp CrankNicolsonDdtScheme::fvcDdtPhiCorr ( @@ -278,6 +277,13 @@ tmp CrankNicolsonDdtScheme::fvcDdtPhiCorr const surfaceScalarField& phi ); +template<> +tmp CrankNicolsonDdtScheme::fvcDdtUfCorr +( + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& Uf +); template<> tmp CrankNicolsonDdtScheme::fvcDdtPhiCorr diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H index ddaae47a87..1e2b840a13 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H @@ -176,7 +176,6 @@ tmp EulerDdtScheme::fvcDdtUfCorr const GeometricField& Uf ); - template<> tmp EulerDdtScheme::fvcDdtPhiCorr ( @@ -184,6 +183,13 @@ tmp EulerDdtScheme::fvcDdtPhiCorr const surfaceScalarField& phi ); +template<> +tmp EulerDdtScheme::fvcDdtUfCorr +( + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& Uf +); template<> tmp EulerDdtScheme::fvcDdtPhiCorr diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H index 8c4f0a178b..b420f2b161 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H @@ -199,7 +199,6 @@ tmp SLTSDdtScheme::fvcDdtUfCorr const GeometricField& Uf ); - template<> tmp SLTSDdtScheme::fvcDdtPhiCorr ( @@ -207,6 +206,13 @@ tmp SLTSDdtScheme::fvcDdtPhiCorr const surfaceScalarField& phi ); +template<> +tmp SLTSDdtScheme::fvcDdtUfCorr +( + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& Uf +); template<> tmp SLTSDdtScheme::fvcDdtPhiCorr diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H index c905f5e29a..425e78ece6 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H @@ -187,7 +187,6 @@ tmp backwardDdtScheme::fvcDdtUfCorr const GeometricField& Uf ); - template<> tmp backwardDdtScheme::fvcDdtPhiCorr ( @@ -195,6 +194,13 @@ tmp backwardDdtScheme::fvcDdtPhiCorr const surfaceScalarField& phi ); +template<> +tmp backwardDdtScheme::fvcDdtUfCorr +( + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& Uf +); template<> tmp backwardDdtScheme::fvcDdtPhiCorr diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H index 2494f8dfcb..db309c8c81 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H @@ -185,7 +185,6 @@ tmp boundedDdtScheme::fvcDdtUfCorr const GeometricField& Uf ); - template<> tmp boundedDdtScheme::fvcDdtPhiCorr ( @@ -193,6 +192,13 @@ tmp boundedDdtScheme::fvcDdtPhiCorr const surfaceScalarField& phi ); +template<> +tmp boundedDdtScheme::fvcDdtUfCorr +( + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& Uf +); template<> tmp boundedDdtScheme::fvcDdtPhiCorr diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H index c78f589122..50aed99d52 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H @@ -91,6 +91,13 @@ public: // Constructors + //- Construct from mesh and name of the rDeltaT field + localEulerDdtScheme(const fvMesh& mesh, const word& rDeltaTName) + : + ddtScheme(mesh), + rDeltaTName_(rDeltaTName) + {} + //- Construct from mesh and Istream localEulerDdtScheme(const fvMesh& mesh, Istream& is) : @@ -188,7 +195,6 @@ tmp localEulerDdtScheme::fvcDdtUfCorr const GeometricField& Uf ); - template<> tmp localEulerDdtScheme::fvcDdtPhiCorr ( @@ -196,6 +202,13 @@ tmp localEulerDdtScheme::fvcDdtPhiCorr const surfaceScalarField& phi ); +template<> +tmp localEulerDdtScheme::fvcDdtUfCorr +( + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& Uf +); template<> tmp localEulerDdtScheme::fvcDdtPhiCorr diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H index ee2d6937a6..4512d05310 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H @@ -175,7 +175,6 @@ tmp steadyStateDdtScheme::fvcDdtUfCorr const GeometricField& Uf ); - template<> tmp steadyStateDdtScheme::fvcDdtPhiCorr ( @@ -183,6 +182,13 @@ tmp steadyStateDdtScheme::fvcDdtPhiCorr const surfaceScalarField& phi ); +template<> +tmp steadyStateDdtScheme::fvcDdtUfCorr +( + const volScalarField& rho, + const volScalarField& U, + const surfaceScalarField& Uf +); template<> tmp steadyStateDdtScheme::fvcDdtPhiCorr From b994d6a11fd686ce21afc08abe41dc8a0b59b2fc Mon Sep 17 00:00:00 2001 From: Henry Date: Sun, 27 Oct 2013 21:00:46 +0000 Subject: [PATCH 4/4] MULES: Added support for predictor-corrector formulation to interFoam and LTSInterFoam --- .../interFoam/LTSInterFoam/LTSInterFoam.C | 5 +- .../interFoam/LTSInterFoam/alphaEqn.H | 39 -- .../interFoam/LTSInterFoam/alphaEqnSubCycle.H | 35 -- .../interFoam/MRFInterFoam/MRFInterFoam.C | 2 +- .../solvers/multiphase/interFoam/alphaEqn.H | 56 +- .../interFoam/interDyMFoam/interDyMFoam.C | 2 +- .../solvers/multiphase/interFoam/interFoam.C | 2 +- .../interMixingFoam/interMixingFoam.C | 2 +- .../porousInterFoam/porousInterFoam.C | 2 +- .../interPhaseChangeFoam/alphaEqn.H | 10 +- .../interPhaseChangeDyMFoam.C | 2 +- .../interPhaseChangeFoam.C | 2 +- src/finiteVolume/Make/files | 1 + .../fvMatrices/solvers/MULES/CMULES.C | 68 +++ .../fvMatrices/solvers/MULES/CMULES.H | 160 ++++++ .../solvers/MULES/CMULESTemplates.C | 493 ++++++++++++++++++ .../fvMatrices/solvers/MULES/IMULES.H | 9 +- .../fvMatrices/solvers/MULES/MULES.C | 19 - .../fvMatrices/solvers/MULES/MULES.H | 65 +-- .../fvMatrices/solvers/MULES/MULESTemplates.C | 430 --------------- .../ras/waterChannel/system/controlDict | 9 +- .../ras/waterChannel/system/fvSchemes | 11 +- .../ras/waterChannel/system/fvSolution | 43 +- .../cavitatingBullet/system/fvSchemes | 1 - 24 files changed, 847 insertions(+), 621 deletions(-) delete mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H delete mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H create mode 100644 src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C create mode 100644 src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H create mode 100644 src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C index 980d1be630..317a18120c 100644 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C @@ -38,7 +38,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" @@ -79,7 +79,10 @@ int main(int argc, char *argv[]) twoPhaseProperties.correct(); + #define LTSSOLVE #include "alphaEqnSubCycle.H" + #undef LTSSOLVE + interface.correct(); turbulence->correct(); diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H deleted file mode 100644 index faae197670..0000000000 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H +++ /dev/null @@ -1,39 +0,0 @@ -{ - word alphaScheme("div(phi,alpha)"); - word alpharScheme("div(phirb,alpha)"); - - surfaceScalarField phic(mag(phi/mesh.magSf())); - phic = min(interface.cAlpha()*phic, max(phic)); - surfaceScalarField phir(phic*interface.nHatf()); - - for (int aCorr=0; aCorr 1) -{ - dimensionedScalar totalDeltaT = runTime.deltaT(); - surfaceScalarField rhoPhiSum - ( - IOobject - ( - "rhoPhiSum", - runTime.timeName(), - mesh - ), - mesh, - dimensionedScalar("0", rhoPhi.dimensions(), 0) - ); - - for - ( - subCycle alphaSubCycle(alpha1, nAlphaSubCycles); - !(++alphaSubCycle).end(); - ) - { - #include "alphaEqn.H" - rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; - } - - rhoPhi = rhoPhiSum; -} -else -{ - #include "alphaEqn.H" -} - -rho == alpha1*rho1 + alpha2*rho2; diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C index 887ef9e34b..0fa1a50ee2 100644 --- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C +++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C @@ -37,7 +37,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H index a2720e20ee..b0dd8ebef2 100644 --- a/applications/solvers/multiphase/interFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interFoam/alphaEqn.H @@ -6,9 +6,39 @@ phic = min(interface.cAlpha()*phic, max(phic)); surfaceScalarField phir(phic*interface.nHatf()); + tmp tphiAlpha; + + if (MULESCorr) + { + fvScalarMatrix alpha1Eqn + ( + #ifdef LTSSOLVE + fv::localEulerDdtScheme(mesh, rDeltaT.name()).fvmDdt(alpha1) + #else + fv::EulerDdtScheme(mesh).fvmDdt(alpha1) + #endif + + fv::gaussConvectionScheme + ( + mesh, + phi, + upwind(mesh, phi) + ).fvmDiv(phi, alpha1) + ); + + alpha1Eqn.solve(); + + Info<< "Phase-1 volume fraction = " + << alpha1.weightedAverage(mesh.Vsc()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; + + tphiAlpha = alpha1Eqn.flux(); + } + for (int aCorr=0; aCorr tphiAlpha0 ( fvc::flux ( @@ -24,12 +54,32 @@ ) ); - MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0); + if (MULESCorr) + { + tphiAlpha0() -= tphiAlpha(); + #ifdef LTSSOLVE + MULES::LTScorrect(alpha1, tphiAlpha0(), 1, 0); + #else + MULES::correct(alpha1, tphiAlpha0(), 1, 0); + #endif + tphiAlpha() += tphiAlpha0(); + } + else + { + tphiAlpha = tphiAlpha0; + + #ifdef LTSSOLVE + MULES::explicitLTSSolve(alpha1, phi, tphiAlpha(), 1, 0); + #else + MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0); + #endif + } alpha2 = 1.0 - alpha1; - rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2; } + rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2; + Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C index 032f2db730..4a673c6a3f 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C @@ -34,7 +34,7 @@ Description #include "fvCFD.H" #include "dynamicFvMesh.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C index fa9b4fffe4..4ce08d27db 100644 --- a/applications/solvers/multiphase/interFoam/interFoam.C +++ b/applications/solvers/multiphase/interFoam/interFoam.C @@ -38,7 +38,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C index 93cb8a0c3d..12bc57d785 100644 --- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C +++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C @@ -31,7 +31,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "threePhaseMixture.H" #include "threePhaseInterfaceProperties.H" diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C index 64938174dc..3423a4654a 100644 --- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C +++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C @@ -39,7 +39,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H index b6e70fc295..45cba01d5b 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H @@ -16,8 +16,14 @@ { fvScalarMatrix alpha1Eqn ( - fvm::ddt(alpha1) - + fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1) + fv::EulerDdtScheme(mesh).fvmDdt(alpha1) + + fv::gaussConvectionScheme + ( + mesh, + phi, + upwind(mesh, phi) + ).fvmDiv(phi, alpha1) + - fvm::Sp(divU, alpha1) == fvm::Sp(vDotvmcAlphal, alpha1) + vDotcAlphal diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C index 630b94534e..31fa619754 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C @@ -43,7 +43,7 @@ Description #include "fvCFD.H" #include "dynamicFvMesh.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "phaseChangeTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C index c1d13ca926..28ef7f0a69 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C @@ -41,7 +41,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "phaseChangeTwoPhaseMixture.H" diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index a99d3e3a31..212a22baf5 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -220,6 +220,7 @@ fields/surfaceFields/surfaceFields.C fvMatrices/fvMatrices.C fvMatrices/fvScalarMatrix/fvScalarMatrix.C fvMatrices/solvers/MULES/MULES.C +fvMatrices/solvers/MULES/CMULES.C fvMatrices/solvers/MULES/IMULES.C fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C new file mode 100644 index 0000000000..a954dccf56 --- /dev/null +++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C @@ -0,0 +1,68 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "CMULES.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void Foam::MULES::correct +( + volScalarField& psi, + surfaceScalarField& phiPsiCorr, + const scalar psiMax, + const scalar psiMin +) +{ + correct + ( + geometricOneField(), + psi, + phiPsiCorr, + zeroField(), zeroField(), + psiMax, psiMin + ); +} + + +void Foam::MULES::LTScorrect +( + volScalarField& psi, + surfaceScalarField& phiPsiCorr, + const scalar psiMax, + const scalar psiMin +) +{ + LTScorrect + ( + geometricOneField(), + psi, + phiPsiCorr, + zeroField(), zeroField(), + psiMax, psiMin + ); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H new file mode 100644 index 0000000000..2786a86aae --- /dev/null +++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H @@ -0,0 +1,160 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + MULES + +Description + CMULES: Multidimensional universal limiter for + explicit corrected implicit solution. + + Solve a convective-only transport equation using an explicit universal + multi-dimensional limiter to correct an implicit conservative bounded + obtained using rigorously bounded schemes such as Euler-implicit in time + upwind in space. + + Parameters are the variable to solve, the normal convective flux and the + actual explicit flux of the variable which is also used to return limited + flux used in the bounded-solution. + +SourceFiles + CMULES.C + CMULESTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef CMULES_H +#define CMULES_H + +#include "MULES.H" +#include "EulerDdtScheme.H" +#include "localEulerDdtScheme.H" +#include "gaussConvectionScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace MULES +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +void correct +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su +); + +template +void correct +( + const RhoType& rho, + volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +); + +void correct +( + volScalarField& psi, + surfaceScalarField& phiCorr, + const scalar psiMax, + const scalar psiMin +); + +template +void LTScorrect +( + const RhoType& rho, + volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +); + +void LTScorrect +( + volScalarField& psi, + surfaceScalarField& phiCorr, + const scalar psiMax, + const scalar psiMin +); + + +template +void limiterCorr +( + scalarField& allLambda, + const RdeltaTType& rDeltaT, + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +); + +template +void limitCorr +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + const volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace MULES +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "CMULESTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C new file mode 100644 index 0000000000..2ff99745da --- /dev/null +++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C @@ -0,0 +1,493 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "CMULES.H" +#include "fvcSurfaceIntegrate.H" +#include "slicedSurfaceFields.H" +#include "wedgeFvPatch.H" +#include "syncTools.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +void Foam::MULES::correct +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su +) +{ + Info<< "MULES: Correcting " << psi.name() << endl; + + const fvMesh& mesh = psi.mesh(); + + scalarField psiIf(psi.size(), 0); + fvc::surfaceIntegrate(psiIf, phiCorr); + + if (mesh.moving()) + { + psi.internalField() = + ( + rho.field()*psi.internalField()*rDeltaT + + Su.field() + - psiIf + )/(rho.field()*rDeltaT - Sp.field()); + } + else + { + psi.internalField() = + ( + rho.field()*psi.internalField()*rDeltaT + + Su.field() + - psiIf + )/(rho.field()*rDeltaT - Sp.field()); + } + + psi.correctBoundaryConditions(); +} + + +template +void Foam::MULES::correct +( + const RhoType& rho, + volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +) +{ + const fvMesh& mesh = psi.mesh(); + const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); + + const dictionary& MULEScontrols = mesh.solverDict(psi.name()); + + label nLimiterIter + ( + readLabel(MULEScontrols.lookup("nLimiterIter")) + ); + + limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter); + correct(rDeltaT, rho, psi, phiCorr, Sp, Su); +} + + +template +void Foam::MULES::LTScorrect +( + const RhoType& rho, + volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +) +{ + const fvMesh& mesh = psi.mesh(); + + const volScalarField& rDeltaT = + mesh.objectRegistry::lookupObject("rSubDeltaT"); + + const dictionary& MULEScontrols = mesh.solverDict(psi.name()); + + label nLimiterIter + ( + readLabel(MULEScontrols.lookup("nLimiterIter")) + ); + + limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter); + correct(rDeltaT, rho, psi, phiCorr, Sp, Su); +} + + +template +void Foam::MULES::limiterCorr +( + scalarField& allLambda, + const RdeltaTType& rDeltaT, + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +) +{ + const scalarField& psiIf = psi; + const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField(); + + const fvMesh& mesh = psi.mesh(); + + const labelUList& owner = mesh.owner(); + const labelUList& neighb = mesh.neighbour(); + tmp tVsc = mesh.Vsc(); + const scalarField& V = tVsc(); + + const scalarField& phiCorrIf = phiCorr; + const surfaceScalarField::GeometricBoundaryField& phiCorrBf = + phiCorr.boundaryField(); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + scalarField& lambdaIf = lambda; + surfaceScalarField::GeometricBoundaryField& lambdaBf = + lambda.boundaryField(); + + scalarField psiMaxn(psiIf.size(), psiMin); + scalarField psiMinn(psiIf.size(), psiMax); + + scalarField sumPhip(psiIf.size(), VSMALL); + scalarField mSumPhim(psiIf.size(), VSMALL); + + forAll(phiCorrIf, facei) + { + label own = owner[facei]; + label nei = neighb[facei]; + + psiMaxn[own] = max(psiMaxn[own], psiIf[nei]); + psiMinn[own] = min(psiMinn[own], psiIf[nei]); + + psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]); + psiMinn[nei] = min(psiMinn[nei], psiIf[own]); + + scalar phiCorrf = phiCorrIf[facei]; + + if (phiCorrf > 0.0) + { + sumPhip[own] += phiCorrf; + mSumPhim[nei] += phiCorrf; + } + else + { + mSumPhim[own] -= phiCorrf; + sumPhip[nei] -= phiCorrf; + } + } + + forAll(phiCorrBf, patchi) + { + const fvPatchScalarField& psiPf = psiBf[patchi]; + const scalarField& phiCorrPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + if (psiPf.coupled()) + { + const scalarField psiPNf(psiPf.patchNeighbourField()); + + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]); + } + } + else + { + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]); + } + } + + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + scalar phiCorrf = phiCorrPf[pFacei]; + + if (phiCorrf > 0.0) + { + sumPhip[pfCelli] += phiCorrf; + } + else + { + mSumPhim[pfCelli] -= phiCorrf; + } + } + } + + psiMaxn = min(psiMaxn, psiMax); + psiMinn = max(psiMinn, psiMin); + + // scalar smooth = 0.5; + // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax); + // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin); + + psiMaxn = + V + *( + (rho.field()*rDeltaT - Sp.field())*psiMaxn + - Su.field() + - rho.field()*psi.internalField()*rDeltaT + ); + + psiMinn = + V + *( + Su.field() + - (rho.field()*rDeltaT - Sp.field())*psiMinn + + rho.field()*psi.internalField()*rDeltaT + ); + + scalarField sumlPhip(psiIf.size()); + scalarField mSumlPhim(psiIf.size()); + + for (int j=0; j 0.0) + { + sumlPhip[own] += lambdaPhiCorrf; + mSumlPhim[nei] += lambdaPhiCorrf; + } + else + { + mSumlPhim[own] -= lambdaPhiCorrf; + sumlPhip[nei] -= lambdaPhiCorrf; + } + } + + forAll(lambdaBf, patchi) + { + scalarField& lambdaPf = lambdaBf[patchi]; + const scalarField& phiCorrfPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + forAll(lambdaPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei]; + + if (lambdaPhiCorrf > 0.0) + { + sumlPhip[pfCelli] += lambdaPhiCorrf; + } + else + { + mSumlPhim[pfCelli] -= lambdaPhiCorrf; + } + } + } + + forAll(sumlPhip, celli) + { + sumlPhip[celli] = + max(min + ( + (sumlPhip[celli] + psiMaxn[celli]) + /(mSumPhim[celli] - SMALL), + 1.0), 0.0 + ); + + mSumlPhim[celli] = + max(min + ( + (mSumlPhim[celli] + psiMinn[celli]) + /(sumPhip[celli] + SMALL), + 1.0), 0.0 + ); + } + + const scalarField& lambdam = sumlPhip; + const scalarField& lambdap = mSumlPhim; + + forAll(lambdaIf, facei) + { + if (phiCorrIf[facei] > 0.0) + { + lambdaIf[facei] = min + ( + lambdaIf[facei], + min(lambdap[owner[facei]], lambdam[neighb[facei]]) + ); + } + else + { + lambdaIf[facei] = min + ( + lambdaIf[facei], + min(lambdam[owner[facei]], lambdap[neighb[facei]]) + ); + } + } + + + forAll(lambdaBf, patchi) + { + fvsPatchScalarField& lambdaPf = lambdaBf[patchi]; + const scalarField& phiCorrfPf = phiCorrBf[patchi]; + const fvPatchScalarField& psiPf = psiBf[patchi]; + + if (isA(mesh.boundary()[patchi])) + { + lambdaPf = 0; + } + else if (psiPf.coupled()) + { + const labelList& pFaceCells = + mesh.boundary()[patchi].faceCells(); + + forAll(lambdaPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + if (phiCorrfPf[pFacei] > 0.0) + { + lambdaPf[pFacei] = + min(lambdaPf[pFacei], lambdap[pfCelli]); + } + else + { + lambdaPf[pFacei] = + min(lambdaPf[pFacei], lambdam[pfCelli]); + } + } + } + else + { + const labelList& pFaceCells = + mesh.boundary()[patchi].faceCells(); + const scalarField& phiCorrPf = phiCorrBf[patchi]; + + forAll(lambdaPf, pFacei) + { + // Limit outlet faces only + if (phiCorrPf[pFacei] > SMALL*SMALL) + { + label pfCelli = pFaceCells[pFacei]; + + if (phiCorrfPf[pFacei] > 0.0) + { + lambdaPf[pFacei] = + min(lambdaPf[pFacei], lambdap[pfCelli]); + } + else + { + lambdaPf[pFacei] = + min(lambdaPf[pFacei], lambdam[pfCelli]); + } + } + } + } + } + + syncTools::syncFaceList(mesh, allLambda, minEqOp()); + } +} + + +template +void Foam::MULES::limitCorr +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + const volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +) +{ + const fvMesh& mesh = psi.mesh(); + + scalarField allLambda(mesh.nFaces(), 1.0); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + limiterCorr + ( + allLambda, + rDeltaT, + rho, + psi, + phiCorr, + Sp, + Su, + psiMax, + psiMin, + nLimiterIter + ); + + phiCorr *= lambda; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H index 1866e8022f..2d6e4a26a0 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H +++ b/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H @@ -25,7 +25,14 @@ Global IMULES Description - IMULES: Multidimensional universal limiter with implicit solution. + IMULES: Multidimensional universal limiter for implicit solution. + + Solve a convective-only transport equation using an explicit universal + multi-dimensional limiter applied to an implicit formulation requiring + iteration to guarantee boundedness. The number of iterations required + to obtain boundedness increases with the Courant number of the simulation. + + It may be more efficient to use CMULES. SourceFiles IMULES.C diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C index f2898b82b8..cb0970b5c4 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C @@ -69,25 +69,6 @@ void Foam::MULES::explicitLTSSolve } -void Foam::MULES::correct -( - volScalarField& psi, - surfaceScalarField& phiPsiCorr, - const scalar psiMax, - const scalar psiMin -) -{ - correct - ( - geometricOneField(), - psi, - phiPsiCorr, - zeroField(), zeroField(), - psiMax, psiMin - ); -} - - void Foam::MULES::limitSum(UPtrList& phiPsiCorrs) { forAll(phiPsiCorrs[0], facei) diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H index b4253f6ee8..f0525ddc9e 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H @@ -25,7 +25,7 @@ Global MULES Description - MULES: Multidimensional universal limiter with explicit solution. + MULES: Multidimensional universal limiter for explicit solution. Solve a convective-only transport equation using an explicit universal multi-dimensional limiter. @@ -126,38 +126,6 @@ void explicitLTSSolve ); -template -void correct -( - const RdeltaTType& rDeltaT, - const RhoType& rho, - volScalarField& psi, - const surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su -); - -template -void correct -( - const RhoType& rho, - volScalarField& psi, - surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin -); - -void correct -( - volScalarField& psi, - surfaceScalarField& phiCorr, - const scalar psiMax, - const scalar psiMin -); - - template void limiter ( @@ -191,41 +159,12 @@ void limit ); -template -void limiterCorr -( - scalarField& allLambda, - const RdeltaTType& rDeltaT, - const RhoType& rho, - const volScalarField& psi, - const surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin, - const label nLimiterIter -); - -template -void limitCorr -( - const RdeltaTType& rDeltaT, - const RhoType& rho, - const volScalarField& psi, - surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin, - const label nLimiterIter -); - - void limitSum(UPtrList& phiPsiCorrs); template void limitSum(SurfaceScalarFieldList& phiPsiCorrs); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace MULES diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C index a0d2e94daa..a7a5282818 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C @@ -138,74 +138,6 @@ void Foam::MULES::explicitLTSSolve } -template -void Foam::MULES::correct -( - const RdeltaTType& rDeltaT, - const RhoType& rho, - volScalarField& psi, - const surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su -) -{ - Info<< "MULES: Correcting " << psi.name() << endl; - - const fvMesh& mesh = psi.mesh(); - - scalarField psiIf(psi.size(), 0); - fvc::surfaceIntegrate(psiIf, phiCorr); - - if (mesh.moving()) - { - psi.internalField() = - ( - rho.field()*psi.internalField()*rDeltaT - + Su.field() - - psiIf - )/(rho.field()*rDeltaT - Sp.field()); - } - else - { - psi.internalField() = - ( - rho.field()*psi.internalField()*rDeltaT - + Su.field() - - psiIf - )/(rho.field()*rDeltaT - Sp.field()); - } - - psi.correctBoundaryConditions(); -} - - -template -void Foam::MULES::correct -( - const RhoType& rho, - volScalarField& psi, - surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin -) -{ - const fvMesh& mesh = psi.mesh(); - const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); - - const dictionary& MULEScontrols = mesh.solverDict(psi.name()); - - label nLimiterIter - ( - readLabel(MULEScontrols.lookup("nLimiterIter")) - ); - - limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter); - correct(rDeltaT, rho, psi, phiCorr, Sp, Su); -} - - template void Foam::MULES::limiter ( @@ -626,368 +558,6 @@ void Foam::MULES::limit } -template -void Foam::MULES::limiterCorr -( - scalarField& allLambda, - const RdeltaTType& rDeltaT, - const RhoType& rho, - const volScalarField& psi, - const surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin, - const label nLimiterIter -) -{ - const scalarField& psiIf = psi; - const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField(); - - const fvMesh& mesh = psi.mesh(); - - const labelUList& owner = mesh.owner(); - const labelUList& neighb = mesh.neighbour(); - tmp tVsc = mesh.Vsc(); - const scalarField& V = tVsc(); - - const scalarField& phiCorrIf = phiCorr; - const surfaceScalarField::GeometricBoundaryField& phiCorrBf = - phiCorr.boundaryField(); - - slicedSurfaceScalarField lambda - ( - IOobject - ( - "lambda", - mesh.time().timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh, - dimless, - allLambda, - false // Use slices for the couples - ); - - scalarField& lambdaIf = lambda; - surfaceScalarField::GeometricBoundaryField& lambdaBf = - lambda.boundaryField(); - - scalarField psiMaxn(psiIf.size(), psiMin); - scalarField psiMinn(psiIf.size(), psiMax); - - scalarField sumPhip(psiIf.size(), VSMALL); - scalarField mSumPhim(psiIf.size(), VSMALL); - - forAll(phiCorrIf, facei) - { - label own = owner[facei]; - label nei = neighb[facei]; - - psiMaxn[own] = max(psiMaxn[own], psiIf[nei]); - psiMinn[own] = min(psiMinn[own], psiIf[nei]); - - psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]); - psiMinn[nei] = min(psiMinn[nei], psiIf[own]); - - scalar phiCorrf = phiCorrIf[facei]; - - if (phiCorrf > 0.0) - { - sumPhip[own] += phiCorrf; - mSumPhim[nei] += phiCorrf; - } - else - { - mSumPhim[own] -= phiCorrf; - sumPhip[nei] -= phiCorrf; - } - } - - forAll(phiCorrBf, patchi) - { - const fvPatchScalarField& psiPf = psiBf[patchi]; - const scalarField& phiCorrPf = phiCorrBf[patchi]; - - const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); - - if (psiPf.coupled()) - { - const scalarField psiPNf(psiPf.patchNeighbourField()); - - forAll(phiCorrPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]); - psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]); - } - } - else - { - forAll(phiCorrPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]); - psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]); - } - } - - forAll(phiCorrPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - scalar phiCorrf = phiCorrPf[pFacei]; - - if (phiCorrf > 0.0) - { - sumPhip[pfCelli] += phiCorrf; - } - else - { - mSumPhim[pfCelli] -= phiCorrf; - } - } - } - - psiMaxn = min(psiMaxn, psiMax); - psiMinn = max(psiMinn, psiMin); - - // scalar smooth = 0.5; - // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax); - // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin); - - psiMaxn = - V - *( - (rho.field()*rDeltaT - Sp.field())*psiMaxn - - Su.field() - - rho.field()*psi.internalField()*rDeltaT - ); - - psiMinn = - V - *( - Su.field() - - (rho.field()*rDeltaT - Sp.field())*psiMinn - + rho.field()*psi.internalField()*rDeltaT - ); - - scalarField sumlPhip(psiIf.size()); - scalarField mSumlPhim(psiIf.size()); - - for (int j=0; j 0.0) - { - sumlPhip[own] += lambdaPhiCorrf; - mSumlPhim[nei] += lambdaPhiCorrf; - } - else - { - mSumlPhim[own] -= lambdaPhiCorrf; - sumlPhip[nei] -= lambdaPhiCorrf; - } - } - - forAll(lambdaBf, patchi) - { - scalarField& lambdaPf = lambdaBf[patchi]; - const scalarField& phiCorrfPf = phiCorrBf[patchi]; - - const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); - - forAll(lambdaPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei]; - - if (lambdaPhiCorrf > 0.0) - { - sumlPhip[pfCelli] += lambdaPhiCorrf; - } - else - { - mSumlPhim[pfCelli] -= lambdaPhiCorrf; - } - } - } - - forAll(sumlPhip, celli) - { - sumlPhip[celli] = - max(min - ( - (sumlPhip[celli] + psiMaxn[celli]) - /(mSumPhim[celli] - SMALL), - 1.0), 0.0 - ); - - mSumlPhim[celli] = - max(min - ( - (mSumlPhim[celli] + psiMinn[celli]) - /(sumPhip[celli] + SMALL), - 1.0), 0.0 - ); - } - - const scalarField& lambdam = sumlPhip; - const scalarField& lambdap = mSumlPhim; - - forAll(lambdaIf, facei) - { - if (phiCorrIf[facei] > 0.0) - { - lambdaIf[facei] = min - ( - lambdaIf[facei], - min(lambdap[owner[facei]], lambdam[neighb[facei]]) - ); - } - else - { - lambdaIf[facei] = min - ( - lambdaIf[facei], - min(lambdam[owner[facei]], lambdap[neighb[facei]]) - ); - } - } - - - forAll(lambdaBf, patchi) - { - fvsPatchScalarField& lambdaPf = lambdaBf[patchi]; - const scalarField& phiCorrfPf = phiCorrBf[patchi]; - const fvPatchScalarField& psiPf = psiBf[patchi]; - - if (isA(mesh.boundary()[patchi])) - { - lambdaPf = 0; - } - else if (psiPf.coupled()) - { - const labelList& pFaceCells = - mesh.boundary()[patchi].faceCells(); - - forAll(lambdaPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - if (phiCorrfPf[pFacei] > 0.0) - { - lambdaPf[pFacei] = - min(lambdaPf[pFacei], lambdap[pfCelli]); - } - else - { - lambdaPf[pFacei] = - min(lambdaPf[pFacei], lambdam[pfCelli]); - } - } - } - else - { - const labelList& pFaceCells = - mesh.boundary()[patchi].faceCells(); - const scalarField& phiCorrPf = phiCorrBf[patchi]; - - forAll(lambdaPf, pFacei) - { - // Limit outlet faces only - if (phiCorrPf[pFacei] > SMALL*SMALL) - { - label pfCelli = pFaceCells[pFacei]; - - if (phiCorrfPf[pFacei] > 0.0) - { - lambdaPf[pFacei] = - min(lambdaPf[pFacei], lambdap[pfCelli]); - } - else - { - lambdaPf[pFacei] = - min(lambdaPf[pFacei], lambdam[pfCelli]); - } - } - } - } - } - - syncTools::syncFaceList(mesh, allLambda, minEqOp()); - } -} - - -template -void Foam::MULES::limitCorr -( - const RdeltaTType& rDeltaT, - const RhoType& rho, - const volScalarField& psi, - surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin, - const label nLimiterIter -) -{ - const fvMesh& mesh = psi.mesh(); - - scalarField allLambda(mesh.nFaces(), 1.0); - - slicedSurfaceScalarField lambda - ( - IOobject - ( - "lambda", - mesh.time().timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh, - dimless, - allLambda, - false // Use slices for the couples - ); - - limiterCorr - ( - allLambda, - rDeltaT, - rho, - psi, - phiCorr, - Sp, - Su, - psiMax, - psiMin, - nLimiterIter - ); - - phiCorr *= lambda; -} - - template void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs) { diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict b/tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict index be315b40ac..05b109778c 100644 --- a/tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict +++ b/tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict @@ -17,7 +17,7 @@ FoamFile application interFoam; -startFrom latestTime; +startFrom startTime; startTime 0; @@ -25,7 +25,7 @@ stopAt endTime; endTime 100; -deltaT 0.01; +deltaT 0.1; writeControl adjustableRunTime; @@ -46,8 +46,9 @@ timePrecision 6; runTimeModifiable yes; adjustTimeStep yes; -maxCo 0.5; -maxAlphaCo 0.5; + +maxCo 6; +maxAlphaCo 6; maxDeltaT 1; functions diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes index 40a4e9640d..11fc30f272 100644 --- a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes +++ b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes @@ -27,12 +27,13 @@ gradSchemes divSchemes { - div(rho*phi,U) Gauss linearUpwind grad(U); - div(phi,alpha) Gauss vanLeer; - div(phirb,alpha) Gauss interfaceCompression; + default none; - div(phi,k) Gauss upwind; - div(phi,omega) $div(phi,k); + div(rho*phi,U) Gauss linearUpwind grad(U); + div(phi,alpha) Gauss vanLeer; + div(phirb,alpha) Gauss interfaceCompression; + + "div\(phi,(k|omega)\)" Gauss upwind; div((muEff*dev(T(grad(U))))) Gauss linear; } diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution index d2428017c4..f6208ef3d2 100644 --- a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution +++ b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution @@ -20,8 +20,16 @@ solvers alpha.water { nAlphaCorr 1; - nAlphaSubCycles 3; + nAlphaSubCycles 1; cAlpha 1; + + MULESCorr yes; + nLimiterIter 3; + + solver PBiCG; + preconditioner DILU; + tolerance 1e-8; + relTol 0; } pcorr @@ -30,9 +38,9 @@ solvers preconditioner { preconditioner GAMG; - tolerance 1e-05; + tolerance 1e-5; relTol 0; - smoother DICGaussSeidel; + smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; @@ -41,31 +49,43 @@ solvers agglomerator faceAreaPair; mergeLevels 1; } - tolerance 1e-10; + tolerance 1e-5; relTol 0; - maxIter 100; + maxIter 50; } p_rgh { - $pcorr; - tolerance 1e-6; - relTol 0.01; + solver GAMG; + tolerance 5e-9; + relTol 0.01; + + smoother GaussSeidel; + nPreSweeps 0; + nPostSweeps 2; + + cacheAgglomeration true; + + nCellsInCoarsestLevel 10; + agglomerator faceAreaPair; + mergeLevels 1; + + maxIter 50; }; p_rghFinal { $p_rgh; - tolerance 1e-6; + tolerance 5e-9; relTol 0; } "(U|k|omega).*" { solver smoothSolver; - smoother GaussSeidel; + smoother symGaussSeidel; nSweeps 1; - tolerance 1e-7; + tolerance 1e-6; relTol 0.1; }; } @@ -84,6 +104,7 @@ relaxationFactors } equations { + ".*" 1; } } diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes index ded24b94b6..bc5571af9e 100644 --- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes +++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes @@ -32,7 +32,6 @@ divSchemes div(phi,k) Gauss linearUpwind grad(k); div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss interfaceCompression; - UD Gauss upwind; div((muEff*dev(T(grad(U))))) Gauss linear; }