From f2cc03bf8dfde5308f8ea2fcad40750385b2eb77 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 22 Mar 2018 11:06:22 +0000 Subject: [PATCH] MULES: Non-uniform limiting and additional form of limit sum MULES and CMULES have been extended so that the limits can be supplied as fields. These arguments are templated so that zeroField, oneField or UniformField can be used in place of a scalar value with no additional overhead. The flux argument has been removed from the unlimited CMULES correct functions in order to make this templating possible. An additional form of limit sum has also been added to MULES. This limits the flux sum by ofsetting in proportion to the phase fraction, rather than by reducing the magnitude of the fluxes with the same sign as the imbalance. The new procedure makes it possible to limit the flux sum in the presence of constraints without encountering a divide by zero. --- .../solvers/multiphase/VoF/alphaEqn.H | 18 +- .../multiphaseMixtureThermo.C | 4 +- .../multiphase/driftFluxFoam/alphaEqn.H | 15 +- .../interFoam/interMixingFoam/alphaEqn.H | 16 +- .../interPhaseChangeFoam/alphaEqn.H | 8 +- .../multiphaseSystem/multiphaseSystem.C | 8 +- .../multiphaseMixture/multiphaseMixture.C | 8 +- .../multiphaseSystem/multiphaseSystem.C | 48 +--- .../twoPhaseSystem/twoPhaseSystem.C | 8 +- .../multiphase/twoLiquidMixingFoam/alphaEqn.H | 10 +- .../twoPhaseSystem/twoPhaseSystem.C | 9 +- .../fields/Fields/uniformField/UniformField.H | 4 +- .../Fields/uniformField/UniformFieldI.H | 70 +++++- src/OpenFOAM/primitives/one/oneI.H | 35 ++- src/OpenFOAM/primitives/zero/zeroI.H | 36 ++- src/finiteVolume/Make/files | 1 - .../fvMatrices/solvers/MULES/CMULES.C | 51 ----- .../fvMatrices/solvers/MULES/CMULES.H | 80 +++++-- .../solvers/MULES/CMULESTemplates.C | 108 +++++++-- .../fvMatrices/solvers/MULES/MULES.C | 58 +++-- .../fvMatrices/solvers/MULES/MULES.H | 107 +++++++-- .../fvMatrices/solvers/MULES/MULESTemplates.C | 212 ++++++++++++++---- 22 files changed, 661 insertions(+), 253 deletions(-) delete mode 100644 src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C diff --git a/applications/solvers/multiphase/VoF/alphaEqn.H b/applications/solvers/multiphase/VoF/alphaEqn.H index ad9ddcb5e5..a9a959fe45 100644 --- a/applications/solvers/multiphase/VoF/alphaEqn.H +++ b/applications/solvers/multiphase/VoF/alphaEqn.H @@ -133,7 +133,15 @@ if (alphaApplyPrevCorr && talphaPhi1Corr0.valid()) { Info<< "Applying the previous iteration compression flux" << endl; - MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0); + MULES::correct + ( + geometricOneField(), + alpha1, + alphaPhi10, + talphaPhi1Corr0.ref(), + oneField(), + zeroField() + ); alphaPhi10 += talphaPhi1Corr0(); } @@ -182,8 +190,8 @@ talphaPhi1Corr.ref(), Sp, (-Sp*alpha1)(), - 1, - 0 + oneField(), + zeroField() ); // Under-relax the correction for all but the 1st corrector @@ -209,8 +217,8 @@ alphaPhi10, Sp, (Su + divU*min(alpha1(), scalar(1)))(), - 1, - 0 + oneField(), + zeroField() ); } diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C index 2263096cff..1eddd3e398 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C @@ -1051,8 +1051,8 @@ void Foam::multiphaseMixtureThermo::solveAlphas alphaPhiCorr, zeroField(), zeroField(), - 1, - 0, + oneField(), + zeroField(), true ); diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H index 3938fde30f..0a64efbcd8 100644 --- a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H +++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H @@ -32,11 +32,12 @@ MULES::correct ( + geometricOneField(), alpha1, alphaPhi, talphaPhiCorr0.ref(), - mixture.alphaMax(), - 0 + UniformField(mixture.alphaMax()), + zeroField() ); alphaPhi += talphaPhiCorr0(); @@ -71,11 +72,12 @@ MULES::correct ( + geometricOneField(), alpha1, talphaPhiUn(), talphaPhiCorr.ref(), - mixture.alphaMax(), - 0 + UniformField(mixture.alphaMax()), + zeroField() ); // Under-relax the correction for all but the 1st corrector @@ -95,11 +97,12 @@ MULES::explicitSolve ( + geometricOneField(), alpha1, phi, alphaPhi, - mixture.alphaMax(), - 0 + UniformField(mixture.alphaMax()), + zeroField() ); } } diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H index bbcd6affd8..e4e0e32f20 100644 --- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H @@ -87,8 +87,8 @@ alphaPhi1, zeroField(), zeroField(), - 1, - 0 + oneField(), + zeroField() ); } else @@ -103,8 +103,8 @@ alphaPhi1, zeroField(), zeroField(), - 1, - 0 + oneField(), + zeroField() ); } @@ -150,8 +150,8 @@ alphaPhi2, zeroField(), zeroField(), - 1, - 0 + oneField(), + zeroField() ); } else @@ -166,8 +166,8 @@ alphaPhi2, zeroField(), zeroField(), - 1, - 0 + oneField(), + zeroField() ); } diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H index 150ccd1f78..e3bbba19d8 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H @@ -78,8 +78,8 @@ divU*(alpha10 - alpha100) - vDotvmcAlphal*alpha10 )(), - 1, - 0 + oneField(), + zeroField() ); // Under-relax the correction for all but the 1st corrector @@ -103,8 +103,8 @@ talphaPhiCorr.ref(), vDotvmcAlphal, (divU*alpha1 + vDotcAlphal)(), - 1, - 0 + oneField(), + zeroField() ); talphaPhi = talphaPhiCorr; diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 7b82d3d1f5..5cc9422af5 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -132,8 +132,8 @@ void Foam::multiphaseSystem::solveAlphas() alphaPhiCorr, zeroField(), zeroField(), - 1, - 0, + oneField(), + zeroField(), true ); @@ -168,9 +168,7 @@ void Foam::multiphaseSystem::solveAlphas() ( geometricOneField(), phase, - alphaPhi, - zeroField(), - zeroField() + alphaPhi ); phase.alphaPhi() = alphaPhi; diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C index 20c8c6e8cf..119e9764c6 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C +++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C @@ -621,8 +621,8 @@ void Foam::multiphaseMixture::solveAlphas alphaPhiCorr, zeroField(), zeroField(), - 1, - 0, + oneField(), + zeroField(), true ); @@ -658,9 +658,7 @@ void Foam::multiphaseMixture::solveAlphas ( geometricOneField(), alpha, - alphaPhi, - zeroField(), - zeroField() + alphaPhi ); rhoPhi_ += alphaPhi*alpha.rho(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 93d7ec3fb5..79aa52e03c 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -69,8 +69,6 @@ void Foam::multiphaseSystem::calcAlphas() void Foam::multiphaseSystem::solveAlphas() { - bool LTS = fv::localEulerDdt::enabled(mesh_); - forAll(phases(), phasei) { phases()[phasei].correctBoundaryConditions(); @@ -138,40 +136,18 @@ void Foam::multiphaseSystem::solveAlphas() phase.correctInflowOutflow(alphaPhiCorr); - if (LTS) - { - MULES::limit - ( - fv::localEulerDdt::localRDeltaT(mesh_), - geometricOneField(), - phase, - phi_, - alphaPhiCorr, - zeroField(), - zeroField(), - phase.alphaMax(), - 0, - true - ); - } - else - { - const scalar rDeltaT = 1.0/mesh_.time().deltaTValue(); - - MULES::limit - ( - rDeltaT, - geometricOneField(), - phase, - phi_, - alphaPhiCorr, - zeroField(), - zeroField(), - phase.alphaMax(), - 0, - true - ); - } + MULES::limit + ( + geometricOneField(), + phase, + phi_, + alphaPhiCorr, + zeroField(), + zeroField(), + UniformField(phase.alphaMax()), + zeroField(), + true + ); } MULES::limitSum(alphaPhiCorrs); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index 05e767253f..796df17ccf 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -274,8 +274,8 @@ void Foam::twoPhaseSystem::solve() alphaPhi10, (alphaSubCycle.index()*Sp)(), (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(), - phase1_.alphaMax(), - 0 + UniformField(phase1_.alphaMax()), + zeroField() ); if (alphaSubCycle.index() == 1) @@ -300,8 +300,8 @@ void Foam::twoPhaseSystem::solve() alphaPhi1, Sp, Su, - phase1_.alphaMax(), - 0 + UniformField(phase1_.alphaMax()), + zeroField() ); phase1_.alphaPhi() = alphaPhi1; diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H index 5344fb75f2..6e9888c1c5 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H @@ -12,7 +12,15 @@ ) ); - MULES::explicitSolve(alpha1, phi, alphaPhi, 1, 0); + MULES::explicitSolve + ( + geometricOneField(), + alpha1, + phi, + alphaPhi, + oneField(), + zeroField() + ); rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index b5ad975720..5e9761b4fe 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -45,6 +45,7 @@ License #include "fixedValueFvsPatchFields.H" #include "blendingMethod.H" #include "HashPtrTable.H" +#include "UniformField.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -464,8 +465,8 @@ void Foam::twoPhaseSystem::solve() alphaPhic10, (alphaSubCycle.index()*Sp)(), (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(), - phase1_.alphaMax(), - 0 + UniformField(phase1_.alphaMax()), + zeroField() ); if (alphaSubCycle.index() == 1) @@ -490,8 +491,8 @@ void Foam::twoPhaseSystem::solve() alphaPhic1, Sp, Su, - phase1_.alphaMax(), - 0 + UniformField(phase1_.alphaMax()), + zeroField() ); phase1_.alphaPhi() = alphaPhic1; diff --git a/src/OpenFOAM/fields/Fields/uniformField/UniformField.H b/src/OpenFOAM/fields/Fields/uniformField/UniformField.H index 140b4b4fd6..3b9f958a2f 100644 --- a/src/OpenFOAM/fields/Fields/uniformField/UniformField.H +++ b/src/OpenFOAM/fields/Fields/uniformField/UniformField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -61,6 +61,8 @@ public: // Member Operators + inline operator Type() const; + inline Type operator[](const label) const; inline UniformField field() const; diff --git a/src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H b/src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H index 9c7012a7ee..9bcf35579a 100644 --- a/src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H +++ b/src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -34,6 +34,13 @@ inline Foam::UniformField::UniformField(const Type& value) {} +template +inline Foam::UniformField::operator Type() const +{ + return value_; +} + + template inline Type Foam::UniformField::operator[](const label) const { @@ -48,4 +55,65 @@ inline Foam::UniformField Foam::UniformField::field() const } +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +inline UniformField min +( + const UniformField& u1, + const UniformField& u2 +) +{ + return UniformField(min(u1.operator Type(), u2.operator Type())); +} + + +template +inline OtherType min(const UniformField& u, const OtherType& o) +{ + return min(u.operator Type(), o); +} + + +template +inline OtherType min(const OtherType& o, const UniformField& u) +{ + return min(o, u.operator Type()); +} + + +template +inline UniformField max +( + const UniformField& u1, + const UniformField& u2 +) +{ + return UniformField(max(u1.operator Type(), u2.operator Type())); +} + + +template +inline OtherType max(const UniformField& u, const OtherType& o) +{ + return max(u.operator Type(), o); +} + + +template +inline OtherType max(const OtherType& o, const UniformField& u) +{ + return max(o, u.operator Type()); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + // ************************************************************************* // diff --git a/src/OpenFOAM/primitives/one/oneI.H b/src/OpenFOAM/primitives/one/oneI.H index 05be7ad120..b85f7caba3 100644 --- a/src/OpenFOAM/primitives/one/oneI.H +++ b/src/OpenFOAM/primitives/one/oneI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -81,6 +81,39 @@ inline const Type& operator/(const Type& t, const one&) return t; } +inline const one& min(const one& o, const one&) +{ + return o; +} + +template +inline Type min(const one&, const Type& t) +{ + return min(scalar(1), t); +} + +template +inline Type min(const Type& t, const one&) +{ + return min(t, scalar(1)); +} + +inline const one& max(const one& o, const one&) +{ + return o; +} + +template +inline Type max(const one&, const Type& t) +{ + return max(scalar(1), t); +} + +template +inline Type max(const Type& t, const one&) +{ + return max(t, scalar(1)); +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/primitives/zero/zeroI.H b/src/OpenFOAM/primitives/zero/zeroI.H index d69255f159..b4bcb0c2be 100644 --- a/src/OpenFOAM/primitives/zero/zeroI.H +++ b/src/OpenFOAM/primitives/zero/zeroI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "zero.H" +#include "scalar.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -89,6 +90,39 @@ inline zero operator/(const zero&, const Type& t) return Zero; } +inline zero min(const zero&, const zero&) +{ + return Zero; +} + +template +inline Type min(const zero&, const Type& t) +{ + return min(scalar(0), t); +} + +template +inline Type min(const Type& t, const zero&) +{ + return min(t, scalar(0)); +} + +inline zero max(const zero&, const zero&) +{ + return Zero; +} + +template +inline Type max(const zero&, const Type& t) +{ + return max(scalar(0), t); +} + +template +inline Type max(const Type& t, const zero&) +{ + return max(t, scalar(0)); +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index da45466621..5f56d4e7f4 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -240,7 +240,6 @@ fields/surfaceFields/surfaceFields.C fvMatrices/fvMatrices.C fvMatrices/fvScalarMatrix/fvScalarMatrix.C fvMatrices/solvers/MULES/MULES.C -fvMatrices/solvers/MULES/CMULES.C fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C interpolation = interpolation/interpolation diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C deleted file mode 100644 index d5cd3fba5d..0000000000 --- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C +++ /dev/null @@ -1,51 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 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, - const surfaceScalarField& phi, - surfaceScalarField& phiPsiCorr, - const scalar psiMax, - const scalar psiMin -) -{ - correct - ( - geometricOneField(), - psi, - phi, - phiPsiCorr, - zeroField(), zeroField(), - psiMax, psiMin - ); -} - - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H index fd76c9f2fd..cffdca571f 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H +++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -66,14 +66,49 @@ void correct const RdeltaTType& rDeltaT, const RhoType& rho, volScalarField& psi, - const surfaceScalarField& phi, const surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su ); +template +void correct +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiCorr +); + template void correct +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su +); + +template +void correct +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiCorr, + const PsiMaxType& psiMax, + const PsiMinType& psiMin +); + +template +< + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> +void correct ( const RhoType& rho, volScalarField& psi, @@ -81,20 +116,19 @@ void correct surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ); -void correct -( - volScalarField& psi, - const surfaceScalarField& phi, - surfaceScalarField& phiCorr, - const scalar psiMax, - const scalar psiMin -); - -template +template +< + class RdeltaTType, + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void limiterCorr ( scalarField& allLambda, @@ -105,11 +139,19 @@ void limiterCorr const surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ); -template +template +< + class RdeltaTType, + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void limitCorr ( const RdeltaTType& rDeltaT, @@ -119,8 +161,8 @@ void limitCorr surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ); diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C index beb0cf66ba..be51f1a9c5 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C @@ -38,7 +38,6 @@ void Foam::MULES::correct const RdeltaTType& rDeltaT, const RhoType& rho, volScalarField& psi, - const surfaceScalarField& phi, const surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su @@ -74,8 +73,67 @@ void Foam::MULES::correct } +template +void Foam::MULES::correct +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiCorr +) +{ + correct(rho, psi, phiCorr, zeroField(), zeroField()); +} + + template void Foam::MULES::correct +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su +) +{ + const fvMesh& mesh = psi.mesh(); + + if (fv::localEulerDdt::enabled(mesh)) + { + const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh); + correct(rDeltaT, rho, psi, phiCorr, Sp, Su); + } + else + { + const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); + correct(rDeltaT, rho, psi, phiCorr, Sp, Su); + } +} + + +template +void Foam::MULES::correct +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiCorr, + const PsiMaxType& psiMax, + const PsiMinType& psiMin +) +{ + correct(rho, psi, phi, phiCorr, zeroField(), zeroField(), psiMax, psiMin); +} + + +template +< + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> +void Foam::MULES::correct ( const RhoType& rho, volScalarField& psi, @@ -83,8 +141,8 @@ void Foam::MULES::correct surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ) { const fvMesh& mesh = psi.mesh(); @@ -105,7 +163,8 @@ void Foam::MULES::correct psiMax, psiMin ); - correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su); + + correct(rDeltaT, rho, psi, phiCorr, Sp, Su); } else { @@ -124,12 +183,20 @@ void Foam::MULES::correct psiMin ); - correct(rDeltaT, rho, psi, phi, phiCorr, Sp, Su); + correct(rDeltaT, rho, psi, phiCorr, Sp, Su); } } -template +template +< + class RdeltaTType, + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void Foam::MULES::limiterCorr ( scalarField& allLambda, @@ -140,8 +207,8 @@ void Foam::MULES::limiterCorr const surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ) { const scalarField& psiIf = psi; @@ -199,8 +266,11 @@ void Foam::MULES::limiterCorr surfaceScalarField::Boundary& lambdaBf = lambda.boundaryFieldRef(); - scalarField psiMaxn(psiIf.size(), psiMin); - scalarField psiMinn(psiIf.size(), psiMax); + scalarField psiMaxn(psiIf.size()); + scalarField psiMinn(psiIf.size()); + + psiMaxn = psiMin; + psiMinn = psiMax; scalarField sumPhip(psiIf.size(), 0.0); scalarField mSumPhim(psiIf.size(), 0.0); @@ -265,8 +335,8 @@ void Foam::MULES::limiterCorr { const label pfCelli = pFaceCells[pFacei]; - psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax); - psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin); + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax[pfCelli]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin[pfCelli]); } } @@ -472,7 +542,15 @@ void Foam::MULES::limiterCorr } -template +template +< + class RdeltaTType, + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void Foam::MULES::limitCorr ( const RdeltaTType& rDeltaT, @@ -482,8 +560,8 @@ void Foam::MULES::limitCorr surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ) { const fvMesh& mesh = psi.mesh(); diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C index 76ba38f969..77ab25c3b1 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C @@ -27,27 +27,6 @@ License // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -void Foam::MULES::explicitSolve -( - volScalarField& psi, - const surfaceScalarField& phi, - surfaceScalarField& phiPsi, - const scalar psiMax, - const scalar psiMin -) -{ - explicitSolve - ( - geometricOneField(), - psi, - phi, - phiPsi, - zeroField(), zeroField(), - psiMax, psiMin - ); -} - - void Foam::MULES::limitSum(UPtrList& phiPsiCorrs) { forAll(phiPsiCorrs[0], facei) @@ -97,4 +76,41 @@ void Foam::MULES::limitSum(UPtrList& phiPsiCorrs) } +void Foam::MULES::limitSum +( + const UPtrList& alphas, + UPtrList& phiPsiCorrs, + const labelHashSet& fixed +) +{ + labelHashSet notFixed(identity(phiPsiCorrs.size())); + notFixed -= fixed; + + forAll(phiPsiCorrs[0], facei) + { + scalar alphaNotFixed = 0, corrNotFixed = 0; + forAllConstIter(labelHashSet, notFixed, iter) + { + alphaNotFixed += alphas[iter.key()][facei]; + corrNotFixed += phiPsiCorrs[iter.key()][facei]; + } + + scalar corrFixed = 0; + forAllConstIter(labelHashSet, fixed, iter) + { + corrFixed += phiPsiCorrs[iter.key()][facei]; + } + + const scalar sumCorr = corrNotFixed + corrFixed; + + const scalar lambda = - sumCorr/alphaNotFixed; + + forAllConstIter(labelHashSet, notFixed, iter) + { + phiPsiCorrs[iter.key()][facei] += lambda*alphas[iter.key()][facei]; + } + } +} + + // ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H index 5ece4fefcb..749fbd2c3c 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -50,6 +50,8 @@ SourceFiles #include "zero.H" #include "zeroField.H" #include "UPtrList.H" +#include "HashSet.H" +#include "UniformField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -71,6 +73,14 @@ void explicitSolve const SuType& Su ); +template +void explicitSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiPsi +); + template void explicitSolve ( @@ -81,7 +91,25 @@ void explicitSolve const SuType& Su ); -template +template +void explicitSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiBD, + surfaceScalarField& phiPsi, + const PsiMaxType& psiMax, + const PsiMinType& psiMin +); + +template +< + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void explicitSolve ( const RhoType& rho, @@ -90,20 +118,19 @@ void explicitSolve surfaceScalarField& phiPsi, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ); -void explicitSolve -( - volScalarField& psi, - const surfaceScalarField& phiBD, - surfaceScalarField& phiPsi, - const scalar psiMax, - const scalar psiMin -); - -template +template +< + class RdeltaTType, + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void limiter ( scalarField& allLambda, @@ -114,11 +141,19 @@ void limiter const surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ); -template +template +< + class RdeltaTType, + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void limit ( const RdeltaTType& rDeltaT, @@ -128,17 +163,51 @@ void limit surfaceScalarField& phiPsi, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin, + const PsiMaxType& psiMax, + const PsiMinType& psiMin, const bool returnCorr ); +template +< + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> +void limit +( + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su, + const PsiMaxType& psiMax, + const PsiMinType& psiMin, + const bool returnCorr +); void limitSum(UPtrList& phiPsiCorrs); template void limitSum(SurfaceScalarFieldList& phiPsiCorrs); +void limitSum +( + const UPtrList& alphas, + UPtrList& phiPsiCorrs, + const labelHashSet& fixed +); + +template +void limitSum +( + const SurfaceScalarFieldList& alphas, + SurfaceScalarFieldList& phiPsiCorrs, + const labelHashSet& fixed +); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C index a7f0c17055..ac4aa598c3 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C @@ -78,6 +78,18 @@ void Foam::MULES::explicitSolve } +template +void Foam::MULES::explicitSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiPsi +) +{ + explicitSolve(rho, psi, phiPsi, zeroField(), zeroField()); +} + + template void Foam::MULES::explicitSolve ( @@ -103,7 +115,39 @@ void Foam::MULES::explicitSolve } -template +template +void Foam::MULES::explicitSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiBD, + surfaceScalarField& phiPsi, + const PsiMaxType& psiMax, + const PsiMinType& psiMin +) +{ + explicitSolve + ( + rho, + psi, + phiBD, + phiPsi, + zeroField(), + zeroField(), + psiMax, + psiMin + ); +} + + +template +< + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void Foam::MULES::explicitSolve ( const RhoType& rho, @@ -112,8 +156,8 @@ void Foam::MULES::explicitSolve surfaceScalarField& phiPsi, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ) { const fvMesh& mesh = psi.mesh(); @@ -123,47 +167,27 @@ void Foam::MULES::explicitSolve if (fv::localEulerDdt::enabled(mesh)) { const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh); - - limit - ( - rDeltaT, - rho, - psi, - phi, - phiPsi, - Sp, - Su, - psiMax, - psiMin, - false - ); - + limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false); explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su); } else { const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); - - limit - ( - rDeltaT, - rho, - psi, - phi, - phiPsi, - Sp, - Su, - psiMax, - psiMin, - false - ); - + limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false); explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su); } } -template +template +< + class RdeltaTType, + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void Foam::MULES::limiter ( scalarField& allLambda, @@ -174,8 +198,8 @@ void Foam::MULES::limiter const surfaceScalarField& phiCorr, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin + const PsiMaxType& psiMax, + const PsiMinType& psiMin ) { const scalarField& psiIf = psi; @@ -235,8 +259,11 @@ void Foam::MULES::limiter scalarField& lambdaIf = lambda; surfaceScalarField::Boundary& lambdaBf = lambda.boundaryFieldRef(); - scalarField psiMaxn(psiIf.size(), psiMin); - scalarField psiMinn(psiIf.size(), psiMax); + scalarField psiMaxn(psiIf.size()); + scalarField psiMinn(psiIf.size()); + + psiMaxn = psiMin; + psiMinn = psiMax; scalarField sumPhiBD(psiIf.size(), 0.0); @@ -307,8 +334,8 @@ void Foam::MULES::limiter { const label pfCelli = pFaceCells[pFacei]; - psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax); - psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin); + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiMax[pfCelli]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiMin[pfCelli]); } } @@ -516,7 +543,15 @@ void Foam::MULES::limiter } -template +template +< + class RdeltaTType, + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> void Foam::MULES::limit ( const RdeltaTType& rDeltaT, @@ -526,8 +561,8 @@ void Foam::MULES::limit surfaceScalarField& phiPsi, const SpType& Sp, const SuType& Su, - const scalar psiMax, - const scalar psiMin, + const PsiMaxType& psiMax, + const PsiMinType& psiMin, const bool returnCorr ) { @@ -595,6 +630,42 @@ void Foam::MULES::limit } +template +< + class RhoType, + class SpType, + class SuType, + class PsiMaxType, + class PsiMinType +> +void Foam::MULES::limit +( + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su, + const PsiMaxType& psiMax, + const PsiMinType& psiMin, + const bool rtnCorr +) +{ + const fvMesh& mesh = psi.mesh(); + + if (fv::localEulerDdt::enabled(mesh)) + { + const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh); + limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, rtnCorr); + } + else + { + const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); + limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, rtnCorr); + } +} + + template void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs) { @@ -631,4 +702,59 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs) } +template +void Foam::MULES::limitSum +( + const SurfaceScalarFieldList& alphas, + SurfaceScalarFieldList& phiPsiCorrs, + const labelHashSet& fixed +) +{ + { + UPtrList alphasInternal(alphas.size()); + forAll(alphas, phasei) + { + alphasInternal.set(phasei, &alphas[phasei]); + } + UPtrList phiPsiCorrsInternal(phiPsiCorrs.size()); + forAll(phiPsiCorrs, phasei) + { + phiPsiCorrsInternal.set(phasei, &phiPsiCorrs[phasei]); + } + + limitSum(alphasInternal, phiPsiCorrsInternal, fixed); + } + + const surfaceScalarField::Boundary& bfld = + phiPsiCorrs[0].boundaryField(); + + forAll(bfld, patchi) + { + if (bfld[patchi].coupled()) + { + UPtrList alphasPatch(alphas.size()); + forAll(alphas, phasei) + { + alphasPatch.set + ( + phasei, + &alphas[phasei].boundaryField()[patchi] + ); + } + UPtrList phiPsiCorrsPatch(phiPsiCorrs.size()); + forAll(phiPsiCorrs, phasei) + { + phiPsiCorrsPatch.set + ( + phasei, + &phiPsiCorrs[phasei].boundaryFieldRef()[patchi] + ); + } + + limitSum(alphasPatch, phiPsiCorrsPatch, fixed); + } + } +} + + // ************************************************************************* //