diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 163955056439cd7615e0acf7e020fb3f026ab643..602affd4e00fbc1644b5decb6faf3f3f3451c448 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/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
index 980d1be6306103ff2aea735fdd2cb49eca276ce3..317a18120c18394070602caed5845f7b186d0999 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/Make/files b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
index 205812e978dd9dfe7ace29f84753167d1186ab39..db8b3af426e048fcb0857ec82650dd8a08868468 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/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
deleted file mode 100644
index faae19767035025223dbd22b6be483d7dc6d6ace..0000000000000000000000000000000000000000
--- 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<nAlphaCorr; aCorr++)
- {
- surfaceScalarField phiAlpha
- (
- fvc::flux
- (
- phi,
- alpha1,
- alphaScheme
- )
- + fvc::flux
- (
- -fvc::flux(-phir, alpha2, alpharScheme),
- alpha1,
- alpharScheme
- )
- );
-
- MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
- //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
-
- alpha2 = 1.0 - alpha1;
- rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
- }
-
- Info<< "Phase-1 volume fraction = "
- << alpha1.weightedAverage(mesh.Vsc()).value()
- << " Min(alpha1) = " << min(alpha1).value()
- << " Max(alpha1) = " << max(alpha1).value()
- << endl;
-}
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
deleted file mode 100644
index 81c65c1caf737480956f72ff2d3952f7311436c5..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
+++ /dev/null
@@ -1,35 +0,0 @@
-#include "alphaControls.H"
-
-if (nAlphaSubCycles > 1)
-{
- dimensionedScalar totalDeltaT = runTime.deltaT();
- surfaceScalarField rhoPhiSum
- (
- IOobject
- (
- "rhoPhiSum",
- runTime.timeName(),
- mesh
- ),
- mesh,
- dimensionedScalar("0", rhoPhi.dimensions(), 0)
- );
-
- for
- (
- subCycle<volScalarField> 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 887ef9e34b5f032d62e730a6fd9faae6a97f0c32..0fa1a50ee2f6613bc8dbd9e4b3b4800b5e0ca9c4 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 a2720e20eeec598b64ac4e783e175670069f98da..b0dd8ebef2abd5bd3e59dc21ce4cc18b8267821c 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<surfaceScalarField> tphiAlpha;
+
+ if (MULESCorr)
+ {
+ fvScalarMatrix alpha1Eqn
+ (
+ #ifdef LTSSOLVE
+ fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
+ #else
+ fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+ #endif
+ + fv::gaussConvectionScheme<scalar>
+ (
+ mesh,
+ phi,
+ upwind<scalar>(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<nAlphaCorr; aCorr++)
{
- surfaceScalarField phiAlpha
+ tmp<surfaceScalarField> 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 032f2db7303cc5004bc0e01355c78bceb346dfa6..4a673c6a3f3770f07b7bfcd5937c4d48254ca9ad 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 fa9b4fffe44c0407fbb51c4639001c272876c75a..4ce08d27db070538ee81555549c40272704a79a4 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/alphaEqns.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H
index 15759adfe4a354816f94c2c928d77558aa9fdd32..970da9fc925c8c15ef0e3aa54f7fde4a2197502f 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/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
index 93cb8a0c3de15eb5e42884b8356aa17019bd1d5f..12bc57d785f7fc5af819bde912225b69c7a3edcf 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 64938174dc06bec138163a2c160a7a176ed2a102..3423a4654affdbe01bea8824f4c693ec0a2b0be9 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 b6e70fc2950b491da1614c4aea58c2e4b80cf955..45cba01d5b62d950cb6b7d795aa1819b3d7502c2 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<scalar>(mesh).fvmDdt(alpha1)
+ + fv::gaussConvectionScheme<scalar>
+ (
+ mesh,
+ phi,
+ upwind<scalar>(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 630b94534eec7a95298de13cda999fa07b8b4ebb..31fa619754e1537253ec4ee1783f8b5d92929730 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 c1d13ca92645a04f75a33c5ce0245e5689461103..28ef7f0a69166534f9860cdd8e55f7ea4fa7c768 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/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 995d07d7e451ac0dd17fd243875fa709ccd2f7eb..0a276c6e6f86a4c8a950f699c30ccddb45657cdb 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 2a5d9a03b0a14c2576ef9485675fab3310eb731d..20afd1f72ffabb375ea915744af2b7c8d0da73f3 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/Make/files b/src/finiteVolume/Make/files
index a99d3e3a3196dd0a5791fdc876c18bde5b85dedd..212a22baf5195e9d20666c4406567f2f40b358d3 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/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C
index 745918403b5528423cc474b6efb6860dfd19ef7c..e5f7a049f4e80181e8033780360a72c7f2801204 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<vector, volMesh>& iF
)
:
- fixedValueFvPatchField<vector>(p, iF),
+ fixedValueFvPatchVectorField(p, iF),
inletVelocity_(p.size(), vector::zero),
alphaName_("alpha")
{}
@@ -53,8 +53,8 @@ interstitialInletVelocityFvPatchVectorField
const fvPatchFieldMapper& mapper
)
:
- fixedValueFvPatchField<vector>(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<vector>(p, iF, dict),
+ fixedValueFvPatchVectorField(p, iF, dict),
inletVelocity_("inletVelocity", dict, p.size()),
alphaName_(dict.lookupOrDefault<word>("alpha", "alpha"))
{}
@@ -79,7 +79,7 @@ interstitialInletVelocityFvPatchVectorField
const interstitialInletVelocityFvPatchVectorField& ptf
)
:
- fixedValueFvPatchField<vector>(ptf),
+ fixedValueFvPatchVectorField(ptf),
inletVelocity_(ptf.inletVelocity_),
alphaName_(ptf.alphaName_)
{}
@@ -92,7 +92,7 @@ interstitialInletVelocityFvPatchVectorField
const DimensionedField<vector, volMesh>& iF
)
:
- fixedValueFvPatchField<vector>(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<const interstitialInletVelocityFvPatchVectorField>(ptf);
+
+ inletVelocity_.rmap(tiptf.inletVelocity_, addr);
+}
+
+
void Foam::interstitialInletVelocityFvPatchVectorField::updateCoeffs()
{
if (updated())
@@ -111,7 +136,7 @@ void Foam::interstitialInletVelocityFvPatchVectorField::updateCoeffs()
patch().lookupPatchField<volScalarField, scalar>(alphaName_);
operator==(inletVelocity_/alphap);
- fixedValueFvPatchField<vector>::updateCoeffs();
+ fixedValueFvPatchVectorField::updateCoeffs();
}
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.H
index 1952b826335ec9556cddb01beab19149f8e71515..6affdefada49a1e8ac53587b62459a1618e4d7e8 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;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
index 58a6d067a94b4bf61b50cd31cf193a51353ded74..204a79358714e9fc2246972e3030e156e2269264 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
@@ -198,7 +198,6 @@ tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
-
template<>
tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr
(
@@ -206,6 +205,13 @@ tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
+template<>
+tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const volScalarField& U,
+ const surfaceScalarField& Uf
+);
template<>
tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
index 41fbbd0c28fc26d409ce9d26a0cae9bab19fb874..94c34b8ea82bd658c54c733f2a819c5352e4774d 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
@@ -270,7 +270,6 @@ tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
-
template<>
tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
(
@@ -278,6 +277,13 @@ tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
+template<>
+tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const volScalarField& U,
+ const surfaceScalarField& Uf
+);
template<>
tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
index ddaae47a872d7187accf9c975aa9153e66de22dd..1e2b840a136f4bc3401d8cc47d280698c85a5e12 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
@@ -176,7 +176,6 @@ tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
-
template<>
tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr
(
@@ -184,6 +183,13 @@ tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
+template<>
+tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const volScalarField& U,
+ const surfaceScalarField& Uf
+);
template<>
tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
index 8c4f0a178b3386d7e9f61bde5d85a8d7562dc765..b420f2b16147f2e590ef9b9c71677c91f37279f1 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
@@ -199,7 +199,6 @@ tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
-
template<>
tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr
(
@@ -207,6 +206,13 @@ tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
+template<>
+tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const volScalarField& U,
+ const surfaceScalarField& Uf
+);
template<>
tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H
index c905f5e29a85928b61c9dda1d63eb804c5868fb3..425e78ece6504d666665e5230b1a58dde7df3c9e 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H
@@ -187,7 +187,6 @@ tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
-
template<>
tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr
(
@@ -195,6 +194,13 @@ tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
+template<>
+tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const volScalarField& U,
+ const surfaceScalarField& Uf
+);
template<>
tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H
index 2494f8dfcb38df690ef731cea619f81e241cfcae..db309c8c8150901e5b83376b68aaabe8d9be4e34 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H
@@ -185,7 +185,6 @@ tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
-
template<>
tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr
(
@@ -193,6 +192,13 @@ tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
+template<>
+tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const volScalarField& U,
+ const surfaceScalarField& Uf
+);
template<>
tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H
index c78f5891228c55b21c433493973a893e0a50b010..50aed99d523ddb8c10751d8c10f43422db1f4fd6 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<Type>(mesh),
+ rDeltaTName_(rDeltaTName)
+ {}
+
//- Construct from mesh and Istream
localEulerDdtScheme(const fvMesh& mesh, Istream& is)
:
@@ -188,7 +195,6 @@ tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
-
template<>
tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
(
@@ -196,6 +202,13 @@ tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
+template<>
+tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const volScalarField& U,
+ const surfaceScalarField& Uf
+);
template<>
tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H
index ee2d6937a6df854cb5ac09fea8198a809c2d2775..4512d0531027df863af9e0ebdc347eca82c4b3c5 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H
@@ -175,7 +175,6 @@ tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
-
template<>
tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr
(
@@ -183,6 +182,13 @@ tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
+template<>
+tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const volScalarField& U,
+ const surfaceScalarField& Uf
+);
template<>
tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C
similarity index 70%
rename from applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
rename to src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C
index 1720a267a0bc252879a8a64d205efb0d3c392bfe..a954dccf561db59269c0abf86ff8973167ffdc9b 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -23,57 +23,42 @@ License
\*---------------------------------------------------------------------------*/
-#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"
+#include "CMULES.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void Foam::MULES::explicitLTSSolve
+void Foam::MULES::correct
(
volScalarField& psi,
- const surfaceScalarField& phi,
- surfaceScalarField& phiPsi,
+ surfaceScalarField& phiPsiCorr,
const scalar psiMax,
const scalar psiMin
)
{
- explicitLTSSolve
+ correct
(
geometricOneField(),
psi,
- phi,
- phiPsi,
+ phiPsiCorr,
zeroField(), zeroField(),
psiMax, psiMin
);
}
-void Foam::MULES::implicitSolve
+void Foam::MULES::LTScorrect
(
volScalarField& psi,
- const surfaceScalarField& phi,
- surfaceScalarField& phiPsi,
+ surfaceScalarField& phiPsiCorr,
const scalar psiMax,
const scalar psiMin
)
{
- implicitSolve
+ LTScorrect
(
geometricOneField(),
psi,
- phi,
- phiPsi,
+ phiPsiCorr,
zeroField(), zeroField(),
psiMax, psiMin
);
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H
similarity index 66%
rename from applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
rename to src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H
index 2c2e282ae2768196ed51c18d9dd0f08da20012dc..2786a86aaebe63b0f3a8cf99ff4d8503cd0edaa7 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.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) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -25,67 +25,77 @@ Global
MULES
Description
- MULES: Multidimensional universal limiter with explicit solution.
+ CMULES: Multidimensional universal limiter for
+ explicit corrected implicit solution.
Solve a convective-only transport equation using an explicit universal
- multi-dimensional limiter.
+ 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
- MULES.C
+ CMULES.C
+ CMULESTemplates.C
\*---------------------------------------------------------------------------*/
-#ifndef MULES_H
-#define MULES_H
+#ifndef CMULES_H
+#define CMULES_H
-#include "volFields.H"
-#include "surfaceFieldsFwd.H"
-#include "primitiveFieldsFwd.H"
-#include "zeroField.H"
-#include "geometricOneField.H"
+#include "MULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "gaussConvectionScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
+namespace MULES
+{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-namespace MULES
-{
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void correct
+(
+ const RdeltaTType& rDeltaT,
+ const RhoType& rho,
+ volScalarField& psi,
+ const surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su
+);
template<class RhoType, class SpType, class SuType>
-void explicitLTSSolve
+void correct
(
const RhoType& rho,
volScalarField& psi,
- const surfaceScalarField& phiBD,
- surfaceScalarField& phiPsi,
+ surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
);
-void explicitLTSSolve
+void correct
(
volScalarField& psi,
- const surfaceScalarField& phiBD,
- surfaceScalarField& phiPsi,
+ surfaceScalarField& phiCorr,
const scalar psiMax,
const scalar psiMin
);
template<class RhoType, class SpType, class SuType>
-void implicitSolve
+void LTScorrect
(
const RhoType& rho,
- volScalarField& gamma,
- const surfaceScalarField& phi,
+ volScalarField& psi,
surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
@@ -93,22 +103,22 @@ void implicitSolve
const scalar psiMin
);
-void implicitSolve
+void LTScorrect
(
- volScalarField& gamma,
- const surfaceScalarField& phi,
+ volScalarField& psi,
surfaceScalarField& phiCorr,
const scalar psiMax,
const scalar psiMin
);
-template<class RhoType, class SpType, class SuType>
-void limiter
+
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void limiterCorr
(
scalarField& allLambda,
+ const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
- const surfaceScalarField& phiBD,
const surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
@@ -117,16 +127,30 @@ void limiter
const label nLimiterIter
);
-} // End namespace MULES
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+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 "MULESTemplates.C"
+# include "CMULESTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
similarity index 55%
rename from applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
rename to src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
index 9185897283a8cf2c0e574ff2583b3b692f0b6946..2ff99745da71a2f463d13b281269de7ce8fef172 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -23,105 +23,49 @@ License
\*---------------------------------------------------------------------------*/
-#include "MULES.H"
-#include "upwind.H"
-#include "uncorrectedSnGrad.H"
-#include "gaussConvectionScheme.H"
-#include "gaussLaplacianScheme.H"
-#include "uncorrectedSnGrad.H"
-#include "surfaceInterpolate.H"
+#include "CMULES.H"
#include "fvcSurfaceIntegrate.H"
#include "slicedSurfaceFields.H"
+#include "wedgeFvPatch.H"
#include "syncTools.H"
-#include "fvCFD.H"
-
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::explicitLTSSolve
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void Foam::MULES::correct
(
+ const RdeltaTType& rDeltaT,
const RhoType& rho,
volScalarField& psi,
- const surfaceScalarField& phi,
- surfaceScalarField& phiPsi,
+ const surfaceScalarField& phiCorr,
const SpType& Sp,
- const SuType& Su,
- const scalar psiMax,
- const scalar psiMin
+ const SuType& Su
)
{
- Info<< "MULES: Solving for " << psi.name() << endl;
+ Info<< "MULES: Correcting " << psi.name() << endl;
const fvMesh& mesh = psi.mesh();
- psi.correctBoundaryConditions();
- surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
-
- surfaceScalarField& phiCorr = phiPsi;
- phiCorr -= phiBD;
-
- scalarField allLambda(mesh.nFaces(), 1.0);
-
- slicedSurfaceScalarField lambda
- (
- IOobject
- (
- "lambda",
- mesh.time().timeName(),
- mesh,
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- mesh,
- dimless,
- allLambda,
- false // Use slices for the couples
- );
-
- limiter
- (
- allLambda,
- rho,
- psi,
- phiBD,
- phiCorr,
- Sp.field(),
- Su.field(),
- psiMax,
- psiMin,
- 3
- );
-
- phiPsi = phiBD + lambda*phiCorr;
-
- scalarField& psiIf = psi;
- const scalarField& psi0 = psi.oldTime();
-
- const volScalarField& rDeltaT =
- mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
-
- psiIf = 0.0;
- fvc::surfaceIntegrate(psiIf, phiPsi);
+ scalarField psiIf(psi.size(), 0);
+ fvc::surfaceIntegrate(psiIf, phiCorr);
if (mesh.moving())
{
- psiIf =
+ psi.internalField() =
(
- mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc()
+ rho.field()*psi.internalField()*rDeltaT
+ Su.field()
- psiIf
- )/(rho*rDeltaT - Sp.field());
+ )/(rho.field()*rDeltaT - Sp.field());
}
else
{
- psiIf =
+ psi.internalField() =
(
- rho.oldTime()*psi0*rDeltaT
+ rho.field()*psi.internalField()*rDeltaT
+ Su.field()
- psiIf
- )/(rho*rDeltaT - Sp.field());
+ )/(rho.field()*rDeltaT - Sp.field());
}
psi.correctBoundaryConditions();
@@ -129,12 +73,11 @@ void Foam::MULES::explicitLTSSolve
template<class RhoType, class SpType, class SuType>
-void Foam::MULES::implicitSolve
+void Foam::MULES::correct
(
const RhoType& rho,
volScalarField& psi,
- const surfaceScalarField& phi,
- surfaceScalarField& phiPsi,
+ surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
@@ -142,176 +85,56 @@ void Foam::MULES::implicitSolve
)
{
const fvMesh& mesh = psi.mesh();
+ const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
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
- );
+ limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+ correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
+}
- CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
- }
- scalarField allLambda(allCoLambda);
- //scalarField allLambda(mesh.nFaces(), 1.0);
+template<class RhoType, class SpType, class SuType>
+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();
- slicedSurfaceScalarField lambda
- (
- IOobject
- (
- "lambda",
- mesh.time().timeName(),
- mesh,
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- mesh,
- dimless,
- allLambda,
- false // Use slices for the couples
- );
+ const volScalarField& rDeltaT =
+ mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
- linear<scalar> CDs(mesh);
- upwind<scalar> UDs(mesh, phi);
- //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
+ const dictionary& MULEScontrols = mesh.solverDict(psi.name());
- fvScalarMatrix psiConvectionDiffusion
+ label nLimiterIter
(
- fvm::ddt(rho, psi)
- + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
- //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
- //.fvmLaplacian(Dpsif, psi)
- - fvm::Sp(Sp, psi)
- - Su
+ readLabel(MULEScontrols.lookup("nLimiterIter"))
);
- surfaceScalarField phiBD(psiConvectionDiffusion.flux());
-
- surfaceScalarField& phiCorr = phiPsi;
- phiCorr -= phiBD;
-
- for (label i=0; i<maxIter; i++)
- {
- if (i != 0 && i < 4)
- {
- allLambda = allCoLambda;
- }
-
- limiter
- (
- allLambda,
- rho,
- psi,
- phiBD,
- phiCorr,
- Sp.field(),
- Su.field(),
- psiMax,
- psiMin,
- nLimiterIter
- );
-
- solve
- (
- psiConvectionDiffusion + fvc::div(lambda*phiCorr),
- MULEScontrols
- );
-
- scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
- scalar minPsi = gMin(psi.internalField());
-
- scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
-
- if (unboundedness < maxUnboundedness)
- {
- break;
- }
- else
- {
- Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
- << " min(" << psi.name() << ") = " << minPsi << endl;
-
- phiBD = psiConvectionDiffusion.flux();
-
- /*
- word gammaScheme("div(phi,gamma)");
- word gammarScheme("div(phirb,gamma)");
-
- const surfaceScalarField& phir =
- mesh.lookupObject<surfaceScalarField>("phir");
-
- phiCorr =
- fvc::flux
- (
- phi,
- psi,
- gammaScheme
- )
- + fvc::flux
- (
- -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
- psi,
- gammarScheme
- )
- - phiBD;
- */
- }
- }
-
- phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
+ limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+ correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
}
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::limiter
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void Foam::MULES::limiterCorr
(
scalarField& allLambda,
+ const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
- const surfaceScalarField& phiBD,
const surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
@@ -323,8 +146,6 @@ void Foam::MULES::limiter
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();
@@ -332,13 +153,6 @@ void Foam::MULES::limiter
tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
const scalarField& V = tVsc();
- const volScalarField& rDeltaT =
- mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
-
- const scalarField& phiBDIf = phiBD;
- const surfaceScalarField::GeometricBoundaryField& phiBDBf =
- phiBD.boundaryField();
-
const scalarField& phiCorrIf = phiCorr;
const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
phiCorr.boundaryField();
@@ -367,8 +181,6 @@ void Foam::MULES::limiter
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);
@@ -383,9 +195,6 @@ void Foam::MULES::limiter
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)
@@ -403,14 +212,13 @@ void Foam::MULES::limiter
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());
+ const scalarField psiPNf(psiPf.patchNeighbourField());
forAll(phiCorrPf, pFacei)
{
@@ -435,8 +243,6 @@ void Foam::MULES::limiter
{
label pfCelli = pFaceCells[pFacei];
- sumPhiBD[pfCelli] += phiBDPf[pFacei];
-
scalar phiCorrf = phiCorrPf[pFacei];
if (phiCorrf > 0.0)
@@ -453,39 +259,30 @@ void Foam::MULES::limiter
psiMaxn = min(psiMaxn, psiMax);
psiMinn = max(psiMinn, psiMin);
- //scalar smooth = 0.5;
- //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
- //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
-
- if (mesh.moving())
- {
- tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
-
- psiMaxn =
- V*((rho*rDeltaT - Sp)*psiMaxn - Su)
- - (V0()*rDeltaT)*rho.oldTime()*psi0
- + sumPhiBD;
+ // scalar smooth = 0.5;
+ // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
+ // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
- 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;
+ psiMaxn =
+ V
+ *(
+ (rho.field()*rDeltaT - Sp.field())*psiMaxn
+ - Su.field()
+ - rho.field()*psi.internalField()*rDeltaT
+ );
- psiMinn =
- V*((rho*rDeltaT)*psi0 - (rho.oldTime()*rDeltaT - Sp)*psiMinn + Su)
- - sumPhiBD;
- }
+ 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<nLimiterIter; j++)
+ for (int j=0; j<nLimiterIter; j++)
{
sumlPhip = 0.0;
mSumlPhim = 0.0;
@@ -533,19 +330,21 @@ void Foam::MULES::limiter
}
}
- forAll (sumlPhip, celli)
+ forAll(sumlPhip, celli)
{
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
);
}
@@ -578,20 +377,57 @@ void Foam::MULES::limiter
{
fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
const scalarField& phiCorrfPf = phiCorrBf[patchi];
+ const fvPatchScalarField& psiPf = psiBf[patchi];
- const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
-
- forAll(lambdaPf, pFacei)
+ if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
{
- label pfCelli = pFaceCells[pFacei];
+ lambdaPf = 0;
+ }
+ else if (psiPf.coupled())
+ {
+ const labelList& pFaceCells =
+ mesh.boundary()[patchi].faceCells();
- if (phiCorrfPf[pFacei] > 0.0)
+ forAll(lambdaPf, pFacei)
{
- lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
+ 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
+ }
+ else
+ {
+ const labelList& pFaceCells =
+ mesh.boundary()[patchi].faceCells();
+ const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+ forAll(lambdaPf, pFacei)
{
- lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
+ // 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]);
+ }
+ }
}
}
}
@@ -601,4 +437,57 @@ void Foam::MULES::limiter
}
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+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 1866e8022f158319487a75ffe26cca638914ea8f..2d6e4a26a0860e196f46fe1e4921373f505d2db3 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/IMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C
index 14749ff9b73d68f9f42f7bfb0b44776f2e63f564..542181e90dd89ac75af8984ed402ba4c43d19444 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 ff5994e670ce942784ddbdb283ebc4247896a603..cb0970b5c4f269b5e1022ea060e2aa27448da675 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
@@ -48,19 +48,21 @@ void Foam::MULES::explicitSolve
}
-void Foam::MULES::correct
+void Foam::MULES::explicitLTSSolve
(
volScalarField& psi,
- surfaceScalarField& phiPsiCorr,
+ const surfaceScalarField& phi,
+ surfaceScalarField& phiPsi,
const scalar psiMax,
const scalar psiMin
)
{
- correct
+ explicitLTSSolve
(
geometricOneField(),
psi,
- phiPsiCorr,
+ phi,
+ phiPsi,
zeroField(), zeroField(),
psiMax, psiMin
);
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
index c1d824b9503a16fce8f27befc168e3a40b754e66..f0525ddc9ec8c73d0cf8d9821f91e1ec29dcbb6b 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.
@@ -60,6 +60,17 @@ namespace MULES
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
+void explicitSolve
+(
+ const RdeltaTType& rDeltaT,
+ const RhoType& rho,
+ volScalarField& psi,
+ const surfaceScalarField& phiPsi,
+ const SpType& Sp,
+ const SuType& Su
+);
+
template<class RhoType, class SpType, class SuType>
void explicitSolve
(
@@ -92,42 +103,34 @@ void explicitSolve
const scalar psiMin
);
-
-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 SpType, class SuType>
-void correct
+void explicitLTSSolve
(
const RhoType& rho,
volScalarField& psi,
- surfaceScalarField& phiCorr,
+ const surfaceScalarField& phiBD,
+ surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
);
-void correct
+void explicitLTSSolve
(
volScalarField& psi,
- surfaceScalarField& phiCorr,
+ const surfaceScalarField& phiBD,
+ surfaceScalarField& phiPsi,
const scalar psiMax,
const scalar psiMin
);
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
void limiter
(
scalarField& allLambda,
+ const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phiBD,
@@ -139,9 +142,10 @@ void limiter
const label nLimiterIter
);
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
void limit
(
+ const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phi,
@@ -155,39 +159,12 @@ void limit
);
-template<class RhoType, class SpType, class SuType>
-void limiterCorr
-(
- scalarField& allLambda,
- 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<class RhoType, class SpType, class SuType>
-void limitCorr
-(
- 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<scalarField>& phiPsiCorrs);
template<class SurfaceScalarFieldList>
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 53b64a5f7b550df925f57ca7ea74406ff81a0602..a7a528281827d5b0e55eb49bd8a77583b93f3983 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -32,9 +32,10 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
+ const RdeltaTType& rDeltaT,
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiPsi,
@@ -48,7 +49,6 @@ void Foam::MULES::explicitSolve
scalarField& psiIf = psi;
const scalarField& psi0 = psi.oldTime();
- const scalar deltaT = mesh.time().deltaTValue();
psiIf = 0.0;
fvc::surfaceIntegrate(psiIf, phiPsi);
@@ -58,19 +58,19 @@ 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();
@@ -82,68 +82,45 @@ void Foam::MULES::explicitSolve
(
const RhoType& rho,
volScalarField& psi,
- const surfaceScalarField& phi,
- surfaceScalarField& phiPsi,
+ const surfaceScalarField& phiPsi,
const SpType& Sp,
- const SuType& Su,
- const scalar psiMax,
- const scalar psiMin
+ const SuType& Su
)
{
- psi.correctBoundaryConditions();
- limit(rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
- explicitSolve(rho, psi, phiPsi, Sp, Su);
+ const fvMesh& mesh = psi.mesh();
+ const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
+ explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}
template<class RhoType, class SpType, class SuType>
-void Foam::MULES::correct
+void Foam::MULES::explicitSolve
(
const RhoType& rho,
volScalarField& psi,
- const surfaceScalarField& phiCorr,
+ const surfaceScalarField& phi,
+ surfaceScalarField& phiPsi,
const SpType& Sp,
- const SuType& Su
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin
)
{
- Info<< "MULES: Correcting " << psi.name() << endl;
-
const fvMesh& mesh = psi.mesh();
-
- const scalar deltaT = mesh.time().deltaTValue();
-
- scalarField psiIf(psi.size(), 0);
- fvc::surfaceIntegrate(psiIf, phiCorr);
-
- if (mesh.moving())
- {
- psi.internalField() =
- (
- rho.field()*psi.internalField()/deltaT
- + Su.field()
- - psiIf
- )/(rho.field()/deltaT - Sp.field());
- }
- else
- {
- psi.internalField() =
- (
- rho.field()*psi.internalField()/deltaT
- + Su.field()
- - psiIf
- )/(rho.field()/deltaT - Sp.field());
- }
-
+ const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
psi.correctBoundaryConditions();
+ limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
+ explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}
template<class RhoType, class SpType, class SuType>
-void Foam::MULES::correct
+void Foam::MULES::explicitLTSSolve
(
const RhoType& rho,
volScalarField& psi,
- surfaceScalarField& phiCorr,
+ const surfaceScalarField& phi,
+ surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
@@ -152,22 +129,20 @@ void Foam::MULES::correct
{
const fvMesh& mesh = psi.mesh();
- const dictionary& MULEScontrols = mesh.solverDict(psi.name());
-
- label nLimiterIter
- (
- readLabel(MULEScontrols.lookup("nLimiterIter"))
- );
+ const volScalarField& rDeltaT =
+ mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
- limitCorr(rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
- correct(rho, psi, phiCorr, Sp, Su);
+ psi.correctBoundaryConditions();
+ limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
+ explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::limiter
(
scalarField& allLambda,
+ const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phiBD,
@@ -190,7 +165,6 @@ void Foam::MULES::limiter
const labelUList& neighb = mesh.neighbour();
tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
const scalarField& V = tVsc();
- const scalar deltaT = mesh.time().deltaTValue();
const scalarField& phiBDIf = phiBD;
const surfaceScalarField::GeometricBoundaryField& phiBDBf =
@@ -321,19 +295,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 +315,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 +325,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 +387,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 +490,10 @@ void Foam::MULES::limiter
}
-template<class RhoType, class SpType, class SuType>
+template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::limit
(
+ const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phi,
@@ -558,6 +535,7 @@ void Foam::MULES::limit
limiter
(
allLambda,
+ rDeltaT,
rho,
psi,
phiBD,
@@ -580,364 +558,6 @@ void Foam::MULES::limit
}
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::limiterCorr
-(
- scalarField& allLambda,
- 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<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
- const scalarField& V = tVsc();
- const scalar deltaT = mesh.time().deltaTValue();
-
- 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()/deltaT - Sp.field())*psiMaxn
- - Su.field()
- - rho.field()*psi.internalField()/deltaT
- );
-
- psiMinn =
- V
- *(
- Su.field()
- - (rho.field()/deltaT - Sp.field())*psiMinn
- + rho.field()*psi.internalField()/deltaT
- );
-
- scalarField sumlPhip(psiIf.size());
- scalarField mSumlPhim(psiIf.size());
-
- for (int j=0; j<nLimiterIter; j++)
- {
- sumlPhip = 0.0;
- mSumlPhim = 0.0;
-
- forAll(lambdaIf, facei)
- {
- label own = owner[facei];
- label nei = neighb[facei];
-
- scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
-
- if (lambdaPhiCorrf > 0.0)
- {
- sumlPhip[own] += lambdaPhiCorrf;
- mSumlPhim[nei] += lambdaPhiCorrf;
- }
- else
- {
- mSumlPhim[own] -= lambdaPhiCorrf;
- sumlPhip[nei] -= lambdaPhiCorrf;
- }
- }
-
- forAll(lambdaBf, patchi)
- {
- scalarField& lambdaPf = lambdaBf[patchi];
- const scalarField& phiCorrfPf = phiCorrBf[patchi];
-
- const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
-
- forAll(lambdaPf, pFacei)
- {
- label pfCelli = pFaceCells[pFacei];
-
- scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
-
- if (lambdaPhiCorrf > 0.0)
- {
- sumlPhip[pfCelli] += lambdaPhiCorrf;
- }
- else
- {
- mSumlPhim[pfCelli] -= lambdaPhiCorrf;
- }
- }
- }
-
- forAll(sumlPhip, celli)
- {
- sumlPhip[celli] =
- max(min
- (
- (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
- 1.0), 0.0
- );
-
- mSumlPhim[celli] =
- max(min
- (
- (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
- 1.0), 0.0
- );
- }
-
- const scalarField& lambdam = sumlPhip;
- const scalarField& lambdap = mSumlPhim;
-
- forAll(lambdaIf, facei)
- {
- if (phiCorrIf[facei] > 0.0)
- {
- lambdaIf[facei] = min
- (
- lambdaIf[facei],
- min(lambdap[owner[facei]], lambdam[neighb[facei]])
- );
- }
- else
- {
- lambdaIf[facei] = min
- (
- lambdaIf[facei],
- min(lambdam[owner[facei]], lambdap[neighb[facei]])
- );
- }
- }
-
-
- forAll(lambdaBf, patchi)
- {
- fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
- const scalarField& phiCorrfPf = phiCorrBf[patchi];
- const fvPatchScalarField& psiPf = psiBf[patchi];
-
- if (isA<wedgeFvPatch>(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<scalar>());
- }
-}
-
-
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::limitCorr
-(
- 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,
- rho,
- psi,
- phiCorr,
- Sp,
- Su,
- psiMax,
- psiMin,
- nLimiterIter
- );
-
- phiCorr *= lambda;
-}
-
-
template<class SurfaceScalarFieldList>
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 be315b40aceac6e99721857ec3b1f36b0585c46d..05b109778c550ad3fa1a6cb1a3a94641829382a8 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 40a4e9640d28f74370feee1d7cd01260e7248968..11fc30f2729251a04dbc6d17b621722ea785fbdb 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 d2428017c40d148cdeba93484b4f9caecbd33c46..f6208ef3d23141f4002c87e5a90a0d4207596194 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 ded24b94b605f04fe5d8e12d4d17c552382b5bf5..bc5571af9ea6e138c703ab98c0f193a4519e17d2 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;
}