zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

MULES: Non-uniform limiting and additional form of limit sum

Browse Source 
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<scalar> 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.
This commit is contained in:
Will Bainbridge
2018-03-22 11:06:22 +00:00
parent ea86dc2cf4
commit f2cc03bf8d
22 changed files with 661 additions and 253 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
applications/solvers/multiphase/VoF/alphaEqn.H
View File

@ -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()
);
}

4
applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
View File

@ -1051,8 +1051,8 @@ void Foam::multiphaseMixtureThermo::solveAlphas
alphaPhiCorr,
zeroField(),
zeroField(),
1,
0,
oneField(),
zeroField(),
true
);

15
applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
View File

@ -32,11 +32,12 @@
MULES::correct
(
geometricOneField(),
alpha1,
alphaPhi,
talphaPhiCorr0.ref(),
mixture.alphaMax(),
0
UniformField<scalar>(mixture.alphaMax()),
zeroField()
);
alphaPhi += talphaPhiCorr0();
@ -71,11 +72,12 @@
MULES::correct
(
geometricOneField(),
alpha1,
talphaPhiUn(),
talphaPhiCorr.ref(),
mixture.alphaMax(),
0
UniformField<scalar>(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<scalar>(mixture.alphaMax()),
zeroField()
);
}
}

16
applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H
View File

@ -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()
);
}

8
applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
View File

@ -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;

8
applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
View File

@ -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;

8
applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
View File

@ -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();

48
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
View File

@ -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<scalar>(phase.alphaMax()),
zeroField(),
true
);
}
MULES::limitSum(alphaPhiCorrs);

8
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
View File

@ -274,8 +274,8 @@ void Foam::twoPhaseSystem::solve()
alphaPhi10,
(alphaSubCycle.index()*Sp)(),
(Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
phase1_.alphaMax(),
0
UniformField<scalar>(phase1_.alphaMax()),
zeroField()
);
if (alphaSubCycle.index() == 1)
@ -300,8 +300,8 @@ void Foam::twoPhaseSystem::solve()
alphaPhi1,
Sp,
Su,
phase1_.alphaMax(),
0
UniformField<scalar>(phase1_.alphaMax()),
zeroField()
);
phase1_.alphaPhi() = alphaPhi1;

10
applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
View File

@ -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;

9
applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
View File

@ -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<scalar>(phase1_.alphaMax()),
zeroField()
);
if (alphaSubCycle.index() == 1)
@ -490,8 +491,8 @@ void Foam::twoPhaseSystem::solve()
alphaPhic1,
Sp,
Su,
phase1_.alphaMax(),
0
UniformField<scalar>(phase1_.alphaMax()),
zeroField()
);
phase1_.alphaPhi() = alphaPhic1;

4
src/OpenFOAM/fields/Fields/uniformField/UniformField.H
View File

@ -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;

70
src/OpenFOAM/fields/Fields/uniformField/UniformFieldI.H
View File

@ -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<Type>::UniformField(const Type& value)
{}
template<class Type>
inline Foam::UniformField<Type>::operator Type() const
{
return value_;
}
template<class Type>
inline Type Foam::UniformField<Type>::operator[](const label) const
{
@ -48,4 +55,65 @@ inline Foam::UniformField<Type> Foam::UniformField<Type>::field() const
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
inline UniformField<Type> min
(
const UniformField<Type>& u1,
const UniformField<Type>& u2
)
{
return UniformField<Type>(min(u1.operator Type(), u2.operator Type()));
}
template<class Type, class OtherType>
inline OtherType min(const UniformField<Type>& u, const OtherType& o)
{
return min(u.operator Type(), o);
}
template<class Type, class OtherType>
inline OtherType min(const OtherType& o, const UniformField<Type>& u)
{
return min(o, u.operator Type());
}
template<class Type>
inline UniformField<Type> max
(
const UniformField<Type>& u1,
const UniformField<Type>& u2
)
{
return UniformField<Type>(max(u1.operator Type(), u2.operator Type()));
}
template<class Type, class OtherType>
inline OtherType max(const UniformField<Type>& u, const OtherType& o)
{
return max(u.operator Type(), o);
}
template<class Type, class OtherType>
inline OtherType max(const OtherType& o, const UniformField<Type>& u)
{
return max(o, u.operator Type());
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

35
src/OpenFOAM/primitives/one/oneI.H
View File

@ -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<class Type>
inline Type min(const one&, const Type& t)
{
return min(scalar(1), t);
}
template<class Type>
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<class Type>
inline Type max(const one&, const Type& t)
{
return max(scalar(1), t);
}
template<class Type>
inline Type max(const Type& t, const one&)
{
return max(t, scalar(1));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

36
src/OpenFOAM/primitives/zero/zeroI.H
View File

@ -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<class Type>
inline Type min(const zero&, const Type& t)
{
return min(scalar(0), t);
}
template<class Type>
inline Type min(const Type& t, const zero&)
{
return min(t, scalar(0));
}
inline zero max(const zero&, const zero&)
{
return Zero;
}
template<class Type>
inline Type max(const zero&, const Type& t)
{
return max(scalar(0), t);
}
template<class Type>
inline Type max(const Type& t, const zero&)
{
return max(t, scalar(0));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

1
src/finiteVolume/Make/files
View File

@ -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

51
src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C
View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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
);
}
// ************************************************************************* //

80
src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H
View File

@ -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<class RhoType>
void correct
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiCorr
);
template<class RhoType, class SpType, class SuType>
void correct
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su
);
template<class RhoType, class PsiMaxType, class PsiMinType>
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<class RdeltaTType, class RhoType, class SpType, class SuType>
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<class RdeltaTType, class RhoType, class SpType, class SuType>
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
);

108
src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
View File

@ -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<class RhoType>
void Foam::MULES::correct
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiCorr
)
{
correct(rho, psi, phiCorr, zeroField(), zeroField());
}
template<class RhoType, class SpType, class SuType>
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<class RhoType, class PsiMaxType, class PsiMinType>
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<class RdeltaTType, class RhoType, class SpType, class SuType>
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<class RdeltaTType, class RhoType, class SpType, class SuType>
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();

58
src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
View File

@ -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<scalarField>& phiPsiCorrs)
{
forAll(phiPsiCorrs[0], facei)
@ -97,4 +76,41 @@ void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
}
void Foam::MULES::limitSum
(
const UPtrList<const scalarField>& alphas,
UPtrList<scalarField>& 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];
}
}
}
// ************************************************************************* //

107
src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
View File

@ -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<class RhoType>
void explicitSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiPsi
);
template<class RhoType, class SpType, class SuType>
void explicitSolve
(
@ -81,7 +91,25 @@ void explicitSolve
const SuType& Su
);
template<class RhoType, class SpType, class SuType>
template<class RhoType, class PsiMaxType, class PsiMinType>
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<class RdeltaTType, class RhoType, class SpType, class SuType>
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<class RdeltaTType, class RhoType, class SpType, class SuType>
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<scalarField>& phiPsiCorrs);
template<class SurfaceScalarFieldList>
void limitSum(SurfaceScalarFieldList& phiPsiCorrs);
void limitSum
(
const UPtrList<const scalarField>& alphas,
UPtrList<scalarField>& phiPsiCorrs,
const labelHashSet& fixed
);
template<class SurfaceScalarFieldList>
void limitSum
(
const SurfaceScalarFieldList& alphas,
SurfaceScalarFieldList& phiPsiCorrs,
const labelHashSet& fixed
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

212
src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
View File

@ -78,6 +78,18 @@ void Foam::MULES::explicitSolve
}
template<class RhoType>
void Foam::MULES::explicitSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiPsi
)
{
explicitSolve(rho, psi, phiPsi, zeroField(), zeroField());
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
@ -103,7 +115,39 @@ void Foam::MULES::explicitSolve
}
template<class RhoType, class SpType, class SuType>
template<class RhoType, class PsiMaxType, class PsiMinType>
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<class RdeltaTType, class RhoType, class SpType, class SuType>
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<class RdeltaTType, class RhoType, class SpType, class SuType>
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<class SurfaceScalarFieldList>
void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
{
@ -631,4 +702,59 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
}
template<class SurfaceScalarFieldList>
void Foam::MULES::limitSum
(
const SurfaceScalarFieldList& alphas,
SurfaceScalarFieldList& phiPsiCorrs,
const labelHashSet& fixed
)
{
{
UPtrList<const scalarField> alphasInternal(alphas.size());
forAll(alphas, phasei)
{
alphasInternal.set(phasei, &alphas[phasei]);
}
UPtrList<scalarField> 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<const scalarField> alphasPatch(alphas.size());
forAll(alphas, phasei)
{
alphasPatch.set
(
phasei,
&alphas[phasei].boundaryField()[patchi]
);
}
UPtrList<scalarField> phiPsiCorrsPatch(phiPsiCorrs.size());
forAll(phiPsiCorrs, phasei)
{
phiPsiCorrsPatch.set
(
phasei,
&phiPsiCorrs[phasei].boundaryFieldRef()[patchi]
);
}
limitSum(alphasPatch, phiPsiCorrsPatch, fixed);
}
}
}
// ************************************************************************* //