diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
index be34ab52ed0638ee53e17931059e32e0b0159389..accc3197ad99e84a3a305077a58302550d7f7f3d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
@@ -1,7 +1,5 @@
{
- label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
- label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+ #include "alphaControls.H"
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
index 1a718434a81f998e90476ddc61d48fc6e24ae6c6..cf509804604443692b63e8a000493d6953ef7f55 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
@@ -1,17 +1,5 @@
#include "readTimeControls.H"
- label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
- label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
-
- if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
- {
- FatalErrorIn(args.executable())
- << "Sub-cycling alpha is only allowed for PISO, "
- "i.e. when the number of outer-correctors = 1"
- << exit(FatalError);
- }
-
bool correctPhi =
pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index a42d9119804c8e61031ee217ec38b8722781f236..ed5faa7c86a515ad69565b24d1269267b1e2ade1 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -56,7 +56,7 @@ int main(int argc, char *argv[])
pimpleControl pimple(mesh);
- #include "readControls.H"
+ #include "readTimeControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "CourantNo.H"
@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
- #include "readControls.H"
+ #include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 73babb08f04ae70d23f2b3ed58831060ec0c771e..13193788219fbe0ba68353f9e4fec7cd1902d5bd 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -85,8 +85,8 @@
if (pimple.finalNonOrthogonalIter())
{
- //p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
- //p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
+ p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+ p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
dgdt =
(
@@ -102,7 +102,7 @@
}
}
- p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+ // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
// Update densities from change in p_rgh
rho1 += psi1*(p_rgh - p_rgh_0);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/readControls.H
deleted file mode 100644
index 87a055d641f44e76ee10348cc43fb17c5d558390..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/compressibleInterFoam/readControls.H
+++ /dev/null
@@ -1,13 +0,0 @@
- #include "readTimeControls.H"
-
- label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
- label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
-
- if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
- {
- FatalErrorIn(args.executable())
- << "Sub-cycling alpha is only allowed for PISO operation, "
- "i.e. when the number of outer-correctors = 1"
- << exit(FatalError);
- }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index cb8b4efc3bbd290c25087d8d39672e0734690d2f..9913595cf296cb67ab33cc3e503ef111cd610475 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,4 +1,2 @@
#include "readTimeControls.H"
-
- int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
- int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));
+ #include "alphaControls.H"
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
index 57c78027a471f712cf06e17486b62fcd05cd3014..428876fd2297f7e39ee228d688bc51307dc49c00 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
index 8525a8bf6738a5b230f524e8ffae6af08063e2b2..f419153506c1ef5962aad3693546bf00ee6a8ae1 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
@@ -130,10 +130,7 @@
<< ", " << gMax(1/rDeltaT.internalField()) << endl;
}
- label nAlphaSubCycles
- (
- readLabel(pimpleDict.lookup("nAlphaSubCycles"))
- );
+ #include "alphaControls.H"
rSubDeltaT = rDeltaT*nAlphaSubCycles;
}
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 57c78027a471f712cf06e17486b62fcd05cd3014..428876fd2297f7e39ee228d688bc51307dc49c00 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H
index 97a09ce017d3b753ddf4a91c7dc37d06e41e46c3..9e3f4f6bd47104f5d2bd52ac09e335fb641a1d37 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H
@@ -1,6 +1,4 @@
-const dictionary& pimpleDict = pimple.dict();
-label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
+#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
index 75015542b554c22e0d6d678a44f3cfe549e1e379..9e5461614be9f75ed7360d7719ef0e5ea71b5e93 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
@@ -4,9 +4,41 @@
surfaceScalarField phir("phir", phic*interface.nHatf());
- for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+ Pair<tmp<volScalarField> > vDotAlphal =
+ twoPhaseProperties->vDotAlphal();
+ const volScalarField& vDotcAlphal = vDotAlphal[0]();
+ const volScalarField& vDotvAlphal = vDotAlphal[1]();
+ const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
+
+ tmp<surfaceScalarField> tphiAlpha;
+
+ if (MULESCorr)
{
- surfaceScalarField phiAlpha
+ fvScalarMatrix alpha1Eqn
+ (
+ fvm::ddt(alpha1)
+ + fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1)
+ ==
+ fvm::Sp(vDotvmcAlphal, alpha1)
+ + vDotcAlphal
+ );
+
+ 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();
+ }
+
+ volScalarField alpha10("alpha10", alpha1);
+
+ for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+ {
+ tmp<surfaceScalarField> tphiAlphaCorr
(
fvc::flux
(
@@ -22,66 +54,52 @@
)
);
- Pair<tmp<volScalarField> > vDotAlphal =
- twoPhaseProperties->vDotAlphal();
- const volScalarField& vDotcAlphal = vDotAlphal[0]();
- const volScalarField& vDotvAlphal = vDotAlphal[1]();
+ if (MULESCorr)
+ {
+ tphiAlphaCorr() -= tphiAlpha();
- volScalarField Sp
- (
- IOobject
+
+ volScalarField alpha100("alpha100", alpha10);
+ alpha10 = alpha1;
+
+ MULES::correct
(
- "Sp",
- runTime.timeName(),
- mesh
- ),
- vDotvAlphal - vDotcAlphal
- );
+ geometricOneField(),
+ alpha1,
+ tphiAlphaCorr(),
+ vDotvmcAlphal,
+ (
+ divU*(alpha10 - alpha100)
+ - vDotvmcAlphal*alpha10
+ )(),
+ 1,
+ 0
+ );
- volScalarField Su
- (
- IOobject
+ tphiAlpha() += tphiAlphaCorr();
+ }
+ else
+ {
+ MULES::explicitSolve
(
- "Su",
- runTime.timeName(),
- mesh
- ),
- // Divergence term is handled explicitly to be
- // consistent with the explicit transport solution
- divU*alpha1
- + vDotcAlphal
- );
+ geometricOneField(),
+ alpha1,
+ phi,
+ tphiAlphaCorr(),
+ vDotvmcAlphal,
+ (divU*alpha1 + vDotcAlphal)(),
+ 1,
+ 0
+ );
- //MULES::explicitSolve
- //(
- // geometricOneField(),
- // alpha1,
- // phi,
- // phiAlpha,
- // Sp,
- // Su,
- // 1,
- // 0
- //);
-
- MULES::implicitSolve
- (
- geometricOneField(),
- alpha1,
- phi,
- phiAlpha,
- Sp,
- Su,
- 1,
- 0
- );
+ tphiAlpha = tphiAlphaCorr;
+ }
alpha2 = 1.0 - alpha1;
- rhoPhi +=
- (runTime.deltaT()/totalDeltaT)
- *(phiAlpha*(rho1 - rho2) + phi*rho2);
}
+ rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
+
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
index 5e7fd27f7528e9e1d38756c4e35b6e5cd52c4897..a4338b907f8c0f633d89919d4896d276d21cc16e 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
@@ -11,21 +11,18 @@ surfaceScalarField rhoPhi
);
{
- const dictionary& pimpleDict = pimple.dict();
-
- label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
-
- label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
+ #include "alphaControls.H"
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
volScalarField divU(fvc::div(phi));
- dimensionedScalar totalDeltaT = runTime.deltaT();
-
if (nAlphaSubCycles > 1)
{
+ dimensionedScalar totalDeltaT = runTime.deltaT();
+ surfaceScalarField rhoPhiSum("rhoPhiSum", rhoPhi);
+
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
@@ -33,7 +30,10 @@ surfaceScalarField rhoPhi
)
{
#include "alphaEqn.H"
+ rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
+
+ rhoPhi = rhoPhiSum;
}
else
{
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index edfc34dfb83d56b9243f8a99b74311b1e435fc56..bc6af055f0c7e39e904101bf8695299e4dbd311b 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -41,7 +41,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
-#include "IMULES.H"
+#include "MULES.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 0cde87042d997622c9745f2edd793dfce867c46a..a19b6a3d9b5f61d9d5fd0028438cf6712f9cd468 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -812,9 +812,8 @@ void Foam::multiphaseSystem::solve()
const Time& runTime = mesh_.time();
- const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
-
- label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
+ const dictionary& alphaControls = mesh_.solverDict(phases_.first().name());
+ label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index b3c4ecb11154126969467151b9862a9beda48bce..76b4dab8dd89cf68c35e4981ba9dc66c882d16c5 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -249,15 +249,12 @@ void Foam::multiphaseMixture::solve()
const Time& runTime = mesh_.time();
- const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
-
- label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
-
- scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha")));
-
-
volScalarField& alpha = phases_.first();
+ const dictionary& alphaControls = mesh_.solverDict(alpha.name());
+ label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+ scalar cAlpha(readScalar(alphaControls.lookup("cAlpha")));
+
if (nAlphaSubCycles > 1)
{
surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
index 57c78027a471f712cf06e17486b62fcd05cd3014..428876fd2297f7e39ee228d688bc51307dc49c00 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index 1abe3ef10d4346899ae7b3875b770374ce543a48..813fdd8f10398b062d20b468e59ce8b7b36ba10d 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,5 +1,4 @@
#include "readTimeControls.H"
+ #include "alphaControls.H"
- int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
- int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));
Switch correctAlpha(pimple.dict().lookup("correctAlpha"));
diff --git a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
index 2e3a56007b22561df858fd9675a9534e928b3c9a..ee9339d0f919689701fc87845e4cc039c3faf253 100644
--- a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
+++ b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -40,6 +40,7 @@ Description
#include "pairPatchAgglomeration.H"
#include "labelListIOList.H"
#include "syncTools.H"
+#include "globalIndex.H"
using namespace Foam;
@@ -53,7 +54,7 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createNamedMesh.H"
- const word dictName("faceAgglomerateDict");
+ const word dictName("viewFactorsDict");
#include "setConstantMeshDictionaryIO.H"
@@ -79,16 +80,16 @@ int main(int argc, char *argv[])
boundary.size()
);
- forAll(boundary, patchId)
- {
- const polyPatch& pp = boundary[patchId];
+ label nCoarseFaces = 0;
- label patchI = pp.index();
- finalAgglom[patchI].setSize(pp.size(), 0);
-
- if (!pp.coupled())
+ forAllConstIter(dictionary, agglomDict, iter)
+ {
+ labelList patchIds = boundary.findIndices(iter().keyword());
+ forAll(patchIds, i)
{
- if (agglomDict.found(pp.name()))
+ label patchI = patchIds[i];
+ const polyPatch& pp = boundary[patchI];
+ if (!pp.coupled())
{
Info << "\nAgglomerating patch : " << pp.name() << endl;
pairPatchAgglomeration agglomObject
@@ -99,12 +100,11 @@ int main(int argc, char *argv[])
agglomObject.agglomerate();
finalAgglom[patchI] =
agglomObject.restrictTopBottomAddressing();
- }
- else
- {
- FatalErrorIn(args.executable())
- << "Patch " << pp.name() << " not found in dictionary: "
- << agglomDict.name() << exit(FatalError);
+
+ if (finalAgglom[patchI].size())
+ {
+ nCoarseFaces += max(finalAgglom[patchI] + 1);
+ }
}
}
}
@@ -144,6 +144,7 @@ int main(int argc, char *argv[])
if (writeAgglom)
{
+ globalIndex index(nCoarseFaces);
volScalarField facesAgglomeration
(
IOobject
@@ -158,14 +159,26 @@ int main(int argc, char *argv[])
dimensionedScalar("facesAgglomeration", dimless, 0)
);
+ label coarsePatchIndex = 0;
forAll(boundary, patchId)
{
- fvPatchScalarField& bFacesAgglomeration =
- facesAgglomeration.boundaryField()[patchId];
-
- forAll(bFacesAgglomeration, j)
+ const polyPatch& pp = boundary[patchId];
+ if (pp.size() > 0)
{
- bFacesAgglomeration[j] = finalAgglom[patchId][j];
+ fvPatchScalarField& bFacesAgglomeration =
+ facesAgglomeration.boundaryField()[patchId];
+
+ forAll(bFacesAgglomeration, j)
+ {
+ bFacesAgglomeration[j] =
+ index.toGlobal
+ (
+ Pstream::myProcNo(),
+ finalAgglom[patchId][j] + coarsePatchIndex
+ );
+ }
+
+ coarsePatchIndex += max(finalAgglom[patchId]) + 1;
}
}
diff --git a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerateDict b/applications/utilities/preProcessing/faceAgglomerate/viewFactorsDict
similarity index 91%
rename from applications/utilities/preProcessing/faceAgglomerate/faceAgglomerateDict
rename to applications/utilities/preProcessing/faceAgglomerate/viewFactorsDict
index 9c164f15092a651ef06fa0b1f871dbee69d2fa82..931511380fd587bdd568d071ae2dfa99535fec0e 100644
--- a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerateDict
+++ b/applications/utilities/preProcessing/faceAgglomerate/viewFactorsDict
@@ -17,6 +17,12 @@ FoamFile
// Write agglomeration as a volScalarField with calculated boundary values
writeFacesAgglomeration true;
+//Debug option
+debug 0;
+
+//Dump connectivity rays
+dumpRays false;
+
// Per patch (wildcard possible) the coarsening level
bottomAir_to_heater
{
diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
index a48fd193f124c43532892acfb27ede445b7a7ac4..ed1d9031ec055d74e0c34168d0a5262474b38539 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
@@ -32,6 +32,19 @@ forAll(patches, patchI)
}
}
+labelList triSurfaceToAgglom(5*nFineFaces);
+
+const triSurface localSurface = triangulate
+(
+ patches,
+ includePatches,
+ finalAgglom,
+ triSurfaceToAgglom,
+ globalNumbering,
+ coarsePatches
+);
+
+
distributedTriSurfaceMesh surfacesMesh
(
IOobject
@@ -43,12 +56,13 @@ distributedTriSurfaceMesh surfacesMesh
IOobject::NO_READ,
IOobject::NO_WRITE
),
- triSurfaceTools::triangulate
- (
- patches,
- includePatches
- ),
+ localSurface,
dict
);
+
+triSurfaceToAgglom.resize(surfacesMesh.size());
+
//surfacesMesh.searchableSurface::write();
+
+surfacesMesh.setField(triSurfaceToAgglom);
diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
index 095436b1a115f3a745ebc540e9d033ff1e172e1f..9cc1ccbcd32d8facff75264f06c8ac4ae714bf3d 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
@@ -2,7 +2,7 @@
// Pre-size by assuming a certain percentage is visible.
// Maximum lenght for dynamicList
-const label maxDynListLength = 10000;
+const label maxDynListLength = 100000;
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
@@ -14,12 +14,17 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
DynamicList<label> startIndex(start.size());
DynamicList<label> endIndex(start.size());
+ DynamicList<label> startAgg(start.size());
+ DynamicList<label> endAgg(start.size());
+
const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
+ const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
const pointField& remoteArea = remoteCoarseSf[procI];
const pointField& remoteFc = remoteCoarseCf[procI];
+ const labelField& remoteAgg = remoteCoarseAgg[procI];
label i = 0;
label j = 0;
@@ -29,6 +34,7 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
const point& fc = myFc[i];
const vector& fA = myArea[i];
+ const label& fAgg = myAgg[i];
for (; j < remoteFc.size(); j++)//
{
@@ -36,26 +42,32 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
const point& remFc = remoteFc[j];
const vector& remA = remoteArea[j];
+ const label& remAgg = remoteAgg[j];
+
const vector& d = remFc - fc;
if (((d & fA) < 0.) && ((d & remA) > 0))
{
- start.append(fc + 0.0001*d);
+ start.append(fc + 0.001*d);
startIndex.append(i);
- end.append(fc + 0.9999*d);
+ startAgg.append(globalNumbering.toGlobal(procI, fAgg));
+ end.append(fc + 0.999*d);
label globalI = globalNumbering.toGlobal(procI, j);
endIndex.append(globalI);
+ endAgg.append(globalNumbering.toGlobal(procI, remAgg));
if (startIndex.size() > maxDynListLength)
{
- break;
+ FatalErrorIn
+ (
+ "shootRays"
+ ) << "Dynamic list need from capacity."
+ << "Actual size maxDynListLength : "
+ << maxDynListLength
+ << abort(FatalError);
}
}
}
}
- if (startIndex.size() > maxDynListLength)
- {
- break;
- }
if (j == remoteFc.size())
{
@@ -63,23 +75,102 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
}
}
- List<pointIndexHit> hitInfo(startIndex.size());
- surfacesMesh.findLine(start, end, hitInfo);
+ }while (returnReduce(i < myFc.size(), orOp<bool>()));
+
+ List<pointIndexHit> hitInfo(startIndex.size());
+ surfacesMesh.findLine(start, end, hitInfo);
+
+ // Return hit global agglo index
+ labelList aggHitIndex;
+ surfacesMesh.getField(hitInfo, aggHitIndex);
+
+ DynamicList<label> dRayIs;
- forAll (hitInfo, rayI)
+ // Collect the rays which has not abstacle in bettween in rayStartFace
+ // and rayEndFace. If the ray hit itself get stored in dRayIs
+ forAll (hitInfo, rayI)
+ {
+ if (!hitInfo[rayI].hit())
+ {
+ rayStartFace.append(startIndex[rayI]);
+ rayEndFace.append(endIndex[rayI]);
+ }
+ else if (aggHitIndex[rayI] == startAgg[rayI])
+ {
+ dRayIs.append(rayI);
+ }
+ }
+
+ start.clear();
+
+
+ // Continue rays which hit themself. If they hit the target
+ // agglomeration are added to rayStartFace and rayEndFace
+
+ bool firstLoop = true;
+ DynamicField<point> startHitItself;
+ DynamicField<point> endHitItself;
+ label iter = 0;
+
+ do
+ {
+ labelField rayIs;
+ rayIs.transfer(dRayIs);
+ dRayIs.clear();
+ forAll (rayIs, rayI)
{
- if (!hitInfo[rayI].hit())
+ const label rayID = rayIs[rayI];
+ label hitIndex = -1;
+
+ if (firstLoop)
{
- rayStartFace.append(startIndex[rayI]);
- rayEndFace.append(endIndex[rayI]);
+ hitIndex = rayIs[rayI];
+ }
+ else
+ {
+ hitIndex = rayI;
+ }
+
+ if (hitInfo[hitIndex].hit())
+ {
+ if (aggHitIndex[hitIndex] == startAgg[rayID])
+ {
+ const vector& endP = end[rayID];
+ const vector& startP = hitInfo[hitIndex].hitPoint();
+ const vector& d = endP - startP;
+
+ startHitItself.append(startP + 0.01*d);
+ endHitItself.append(startP + 1.01*d);
+
+ dRayIs.append(rayID);
+ }
+ else if (aggHitIndex[hitIndex] == endAgg[rayID])
+ {
+ rayStartFace.append(startIndex[rayID]);
+ rayEndFace.append(endIndex[rayID]);
+ }
+
}
}
- start.clear();
- startIndex.clear();
- end.clear();
- endIndex.clear();
+ hitInfo.clear();
+ hitInfo.resize(dRayIs.size());
- }while (returnReduce(i < myFc.size(), orOp<bool>()));
+ surfacesMesh.findLine(startHitItself, endHitItself, hitInfo);
+
+ surfacesMesh.getField(hitInfo, aggHitIndex);
+
+
+ endHitItself.clear();
+ startHitItself.clear();
+ firstLoop = false;
+ iter ++;
+
+ }while (returnReduce(hitInfo.size(), orOp<bool>()) > 0 && iter < 10);
+ startIndex.clear();
+ end.clear();
+ endIndex.clear();
+ startAgg.clear();
+ endAgg.clear();
}
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 140ac3da8aaa19a9d293c18e0600e01b143a9497..c9aa219d13227b7b5d4eb56c7ac1e53964017ded 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -68,6 +68,101 @@ Description
using namespace Foam;
+
+triSurface triangulate
+(
+ const polyBoundaryMesh& bMesh,
+ const labelHashSet& includePatches,
+ const labelListIOList& finalAgglom,
+ labelList& triSurfaceToAgglom,
+ const globalIndex& globalNumbering,
+ const polyBoundaryMesh& coarsePatches
+)
+{
+ const polyMesh& mesh = bMesh.mesh();
+
+ // Storage for surfaceMesh. Size estimate.
+ DynamicList<labelledTri> triangles
+ (
+ mesh.nFaces() - mesh.nInternalFaces()
+ );
+
+ label newPatchI = 0;
+ label localTriFaceI = 0;
+
+ forAllConstIter(labelHashSet, includePatches, iter)
+ {
+ const label patchI = iter.key();
+ const polyPatch& patch = bMesh[patchI];
+ const pointField& points = patch.points();
+
+ label nTriTotal = 0;
+
+ forAll(patch, patchFaceI)
+ {
+ const face& f = patch[patchFaceI];
+
+ faceList triFaces(f.nTriangles(points));
+
+ label nTri = 0;
+
+ f.triangles(points, nTri, triFaces);
+
+ forAll(triFaces, triFaceI)
+ {
+ const face& f = triFaces[triFaceI];
+
+ triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
+
+ nTriTotal++;
+
+ triSurfaceToAgglom[localTriFaceI++] = globalNumbering.toGlobal
+ (
+ Pstream::myProcNo(),
+ finalAgglom[patchI][patchFaceI]
+ + coarsePatches[patchI].start()
+ );
+ }
+ }
+
+ newPatchI++;
+ }
+
+ triSurfaceToAgglom.resize(localTriFaceI-1);
+
+ triangles.shrink();
+
+ // Create globally numbered tri surface
+ triSurface rawSurface(triangles, mesh.points());
+
+ // Create locally numbered tri surface
+ triSurface surface
+ (
+ rawSurface.localFaces(),
+ rawSurface.localPoints()
+ );
+
+ // Add patch names to surface
+ surface.patches().setSize(newPatchI);
+
+ newPatchI = 0;
+
+ forAllConstIter(labelHashSet, includePatches, iter)
+ {
+ const label patchI = iter.key();
+ const polyPatch& patch = bMesh[patchI];
+
+ surface.patches()[newPatchI].index() = patchI;
+ surface.patches()[newPatchI].name() = patch.name();
+ surface.patches()[newPatchI].geometricType() = patch.type();
+
+ newPatchI++;
+ }
+
+ return surface;
+}
+
+
void writeRays
(
const fileName& fName,
@@ -213,7 +308,7 @@ int main(int argc, char *argv[])
if (debug)
{
- Info << "\nCreating single cell mesh..." << endl;
+ Pout << "\nCreating single cell mesh..." << endl;
}
singleCellFvMesh coarseMesh
@@ -230,6 +325,11 @@ int main(int argc, char *argv[])
finalAgglom
);
+ if (debug)
+ {
+ Pout << "\nCreated single cell mesh..." << endl;
+ }
+
// Calculate total number of fine and coarse faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -283,6 +383,8 @@ int main(int argc, char *argv[])
DynamicList<point> localCoarseCf(nCoarseFaces);
DynamicList<point> localCoarseSf(nCoarseFaces);
+ DynamicList<label> localAgg(nCoarseFaces);
+
forAll (viewFactorsPatches, i)
{
const label patchID = viewFactorsPatches[i];
@@ -296,11 +398,18 @@ int main(int argc, char *argv[])
const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchID];
const pointField& coarseSf = coarseMesh.Sf().boundaryField()[patchID];
+ labelHashSet includePatches;
+ includePatches.insert(patchID);
+
forAll(coarseCf, faceI)
{
point cf = coarseCf[faceI];
+
const label coarseFaceI = coarsePatchFace[faceI];
const labelList& fineFaces = coarseToFine[coarseFaceI];
+ const label agglomI =
+ agglom[fineFaces[0]] + coarsePatches[patchID].start();
+
// Construct single face
uindirectPrimitivePatch upp
(
@@ -308,6 +417,7 @@ int main(int argc, char *argv[])
pp.points()
);
+
List<point> availablePoints
(
upp.faceCentres().size()
@@ -342,6 +452,7 @@ int main(int argc, char *argv[])
point sf = coarseSf[faceI];
localCoarseCf.append(cf);
localCoarseSf.append(sf);
+ localAgg.append(agglomI);
}
}
@@ -350,9 +461,12 @@ int main(int argc, char *argv[])
List<pointField> remoteCoarseCf(Pstream::nProcs());
List<pointField> remoteCoarseSf(Pstream::nProcs());
+ List<labelField> remoteCoarseAgg(Pstream::nProcs());
remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
+ remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
+
// Collect remote Cf and Sf on fine mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -370,8 +484,12 @@ int main(int argc, char *argv[])
Pstream::scatterList(remoteCoarseCf);
Pstream::gatherList(remoteCoarseSf);
Pstream::scatterList(remoteCoarseSf);
+ Pstream::gatherList(remoteCoarseAgg);
+ Pstream::scatterList(remoteCoarseAgg);
+ globalIndex globalNumbering(nCoarseFaces);
+
// Set up searching engine for obstacles
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#include "searchingEngine.H"
@@ -383,8 +501,6 @@ int main(int argc, char *argv[])
DynamicList<label> rayEndFace(rayStartFace.size());
- globalIndex globalNumbering(nCoarseFaces);
-
// Return rayStartFace in local index andrayEndFace in global index
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/applications/utilities/surface/surfaceHookUp/Make/files b/applications/utilities/surface/surfaceHookUp/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..b495b65457b4eee142d4de5cdaf1b1ee607c87e6
--- /dev/null
+++ b/applications/utilities/surface/surfaceHookUp/Make/files
@@ -0,0 +1,3 @@
+surfaceHookUp.C
+
+EXE = $(FOAM_APPBIN)/surfaceHookUp
diff --git a/applications/utilities/surface/surfaceHookUp/Make/options b/applications/utilities/surface/surfaceHookUp/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..db53089d48087fcbd6df6c333a0527abbfea4609
--- /dev/null
+++ b/applications/utilities/surface/surfaceHookUp/Make/options
@@ -0,0 +1,9 @@
+EXE_INC = \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+ -I$(LIB_SRC)/triSurface/lnInclude
+
+EXE_LIBS = \
+ -lmeshTools \
+ -lfileFormats \
+ -ltriSurface
diff --git a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C
new file mode 100644
index 0000000000000000000000000000000000000000..61829b404b971f4e3dd2211b201ced2a981bfc45
--- /dev/null
+++ b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C
@@ -0,0 +1,525 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Application
+ surfaceHookUp
+
+Description
+ Find close open edges and stitches the surface along them
+
+Usage
+ - surfaceHookUp hookDistance [OPTION]
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+
+#include "triSurfaceMesh.H"
+#include "indexedOctree.H"
+#include "treeBoundBox.H"
+#include "PackedBoolList.H"
+#include "unitConversion.H"
+#include "searchableSurfaces.H"
+
+using namespace Foam;
+
+// Split faceI along edgeI at position newPointI
+void greenRefine
+(
+ const triSurface& surf,
+ const label faceI,
+ const label edgeI,
+ const label newPointI,
+ DynamicList<labelledTri>& newFaces
+)
+{
+ const labelledTri& f = surf.localFaces()[faceI];
+ const edge& e = surf.edges()[edgeI];
+
+ // Find index of edge in face.
+
+ label fp0 = findIndex(f, e[0]);
+ label fp1 = f.fcIndex(fp0);
+ label fp2 = f.fcIndex(fp1);
+
+ if (f[fp1] == e[1])
+ {
+ // Edge oriented like face
+ newFaces.append
+ (
+ labelledTri
+ (
+ f[fp0],
+ newPointI,
+ f[fp2],
+ f.region()
+ )
+ );
+ newFaces.append
+ (
+ labelledTri
+ (
+ newPointI,
+ f[fp1],
+ f[fp2],
+ f.region()
+ )
+ );
+ }
+ else
+ {
+ newFaces.append
+ (
+ labelledTri
+ (
+ f[fp2],
+ newPointI,
+ f[fp1],
+ f.region()
+ )
+ );
+ newFaces.append
+ (
+ labelledTri
+ (
+ newPointI,
+ f[fp0],
+ f[fp1],
+ f.region()
+ )
+ );
+ }
+}
+
+
+//scalar checkEdgeAngle
+//(
+// const triSurface& surf,
+// const label edgeIndex,
+// const label pointIndex,
+// const scalar& angle
+//)
+//{
+// const edge& e = surf.edges()[edgeIndex];
+
+// vector eVec = e.vec(surf.localPoints());
+// eVec /= mag(eVec) + SMALL;
+
+// const labelList& pEdges = surf.pointEdges()[pointIndex];
+//
+// forAll(pEdges, eI)
+// {
+// const edge& nearE = surf.edges()[pEdges[eI]];
+
+// vector nearEVec = nearE.vec(surf.localPoints());
+// nearEVec /= mag(nearEVec) + SMALL;
+
+// const scalar dot = eVec & nearEVec;
+// const scalar minCos = degToRad(angle);
+
+// if (mag(dot) > minCos)
+// {
+// return false;
+// }
+// }
+
+// return true;
+//}
+
+
+void createBoundaryEdgeTrees
+(
+ const PtrList<triSurfaceMesh>& surfs,
+ PtrList<indexedOctree<treeDataEdge> >& bEdgeTrees,
+ labelListList& treeBoundaryEdges
+)
+{
+ forAll(surfs, surfI)
+ {
+ const triSurface& surf = surfs[surfI];
+
+ // Boundary edges
+ treeBoundaryEdges[surfI] =
+ labelList
+ (
+ identity(surf.nEdges() - surf.nInternalEdges())
+ + surf.nInternalEdges()
+ );
+
+ Random rndGen(17301893);
+
+ // Slightly extended bb. Slightly off-centred just so on symmetric
+ // geometry there are less face/edge aligned items.
+ treeBoundBox bb
+ (
+ treeBoundBox(UList<point>(surf.localPoints())).extend(rndGen, 1e-4)
+ );
+
+ bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+ bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+ bEdgeTrees.set
+ (
+ surfI,
+ new indexedOctree<treeDataEdge>
+ (
+ treeDataEdge
+ (
+ false, // cachebb
+ surf.edges(), // edges
+ surf.localPoints(), // points
+ treeBoundaryEdges[surfI] // selected edges
+ ),
+ bb, // bb
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
+ )
+ );
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ argList::addNote
+ (
+ "hook surfaces to other surfaces by moving and retriangulating their"
+ "boundary edges to match other surface boundary edges"
+ );
+ argList::noParallel();
+ argList::validArgs.append("hookTolerance");
+
+# include "addDictOption.H"
+
+# include "setRootCase.H"
+# include "createTime.H"
+
+ const word dictName("surfaceHookUpDict");
+# include "setSystemRunTimeDictionaryIO.H"
+
+ Info<< "Reading " << dictName << nl << endl;
+
+ const IOdictionary dict(dictIO);
+
+ const scalar dist(args.argRead<scalar>(1));
+ const scalar matchTolerance(SMALL);
+
+ Info<< "Hooking distance = " << dist << endl;
+
+ searchableSurfaces surfs
+ (
+ IOobject
+ (
+ "surfacesToHook",
+ runTime.constant(),
+ "triSurface",
+ runTime
+ ),
+ dict
+ );
+
+ Info<< nl << "Reading surfaces: " << endl;
+ forAll(surfs, surfI)
+ {
+ Info<< incrIndent;
+ Info<< nl << indent << "Surface = " << surfs.names()[surfI] << endl;
+
+ const wordList& regions = surfs[surfI].regions();
+ forAll(regions, surfRegionI)
+ {
+ Info<< incrIndent;
+ Info<< indent << "Regions = " << regions[surfRegionI] << endl;
+ Info<< decrIndent;
+ }
+ Info<< decrIndent;
+ }
+
+ PtrList<indexedOctree<treeDataEdge> > bEdgeTrees(surfs.size());
+ labelListList treeBoundaryEdges(surfs.size());
+
+ List<DynamicList<labelledTri> > newFaces(surfs.size());
+ List<DynamicList<point> > newPoints(surfs.size());
+ List<PackedBoolList> visitedFace(surfs.size());
+
+ PtrList<triSurfaceMesh> newSurfaces(surfs.size());
+ forAll(surfs, surfI)
+ {
+ const triSurfaceMesh& surf =
+ refCast<const triSurfaceMesh>(surfs[surfI]);
+
+ newSurfaces.set
+ (
+ surfI,
+ new triSurfaceMesh
+ (
+ IOobject
+ (
+ "hookedSurface_" + surfs.names()[surfI],
+ runTime.constant(),
+ "triSurface",
+ runTime
+ ),
+ surf
+ )
+ );
+ }
+
+ label nChanged = 0;
+ label nIters = 0;
+
+ do
+ {
+ Info<< nl << "Iteration = " << nIters++ << endl;
+ nChanged = 0;
+
+ createBoundaryEdgeTrees(newSurfaces, bEdgeTrees, treeBoundaryEdges);
+
+ forAll(newSurfaces, surfI)
+ {
+ const triSurface& newSurf = newSurfaces[surfI];
+
+ newFaces[surfI] = newSurf.localFaces();
+ newPoints[surfI] = newSurf.localPoints();
+ visitedFace[surfI] = PackedBoolList(newSurf.size(), false);
+ }
+
+ forAll(newSurfaces, surfI)
+ {
+ const triSurface& surf = newSurfaces[surfI];
+
+ List<pointIndexHit> bPointsTobEdges(surf.boundaryPoints().size());
+ labelList bPointsHitTree(surf.boundaryPoints().size(), -1);
+
+ const labelListList& pointEdges = surf.pointEdges();
+
+ forAll(bPointsTobEdges, bPointI)
+ {
+ pointIndexHit& nearestHit = bPointsTobEdges[bPointI];
+
+ const label pointI = surf.boundaryPoints()[bPointI];
+ const point& samplePt = surf.localPoints()[pointI];
+
+ const labelList& pEdges = pointEdges[pointI];
+
+ // Add edges connected to the edge to the shapeMask
+ DynamicList<label> shapeMask;
+ shapeMask.append(pEdges);
+
+ forAll(bEdgeTrees, treeI)
+ {
+ const indexedOctree<treeDataEdge>& bEdgeTree =
+ bEdgeTrees[treeI];
+
+ pointIndexHit currentHit =
+ bEdgeTree.findNearest
+ (
+ samplePt,
+ sqr(dist),
+ treeDataEdge::findNearestOpSubset
+ (
+ bEdgeTree,
+ shapeMask
+ )
+ );
+
+ if
+ (
+ currentHit.hit()
+ &&
+ (
+ !nearestHit.hit()
+ ||
+ (
+ magSqr(currentHit.hitPoint() - samplePt)
+ < magSqr(nearestHit.hitPoint() - samplePt)
+ )
+ )
+ )
+ {
+ nearestHit = currentHit;
+ bPointsHitTree[bPointI] = treeI;
+ }
+ }
+
+ scalar dist2 = magSqr(nearestHit.rawPoint() - samplePt);
+
+ if (nearestHit.hit())
+ {
+ // bool rejectEdge =
+ // checkEdgeAngle
+ // (
+ // surf,
+ // nearestHit.index(),
+ // pointI,
+ // 30
+ // );
+
+ if (dist2 > Foam::sqr(dist))
+ {
+ nearestHit.setMiss();
+ }
+ }
+ }
+
+ forAll(bPointsTobEdges, bPointI)
+ {
+ const pointIndexHit& eHit = bPointsTobEdges[bPointI];
+
+ if (eHit.hit())
+ {
+ const label hitSurfI = bPointsHitTree[bPointI];
+ const triSurface& hitSurf = newSurfaces[hitSurfI];
+
+ const label eIndex =
+ treeBoundaryEdges[hitSurfI][eHit.index()];
+ const edge& e = hitSurf.edges()[eIndex];
+
+ const label pointI = surf.boundaryPoints()[bPointI];
+
+ const labelList& eFaces = hitSurf.edgeFaces()[eIndex];
+
+ if (eFaces.size() != 1)
+ {
+ WarningIn(args.executable())
+ << "Edge is attached to " << eFaces.size()
+ << " faces." << endl;
+
+ continue;
+ }
+
+ const label faceI = eFaces[0];
+
+ if (visitedFace[hitSurfI][faceI])
+ {
+ continue;
+ }
+
+ DynamicList<labelledTri> newFacesFromSplit(2);
+
+ const point& pt = surf.localPoints()[pointI];
+
+ if
+ (
+ (
+ magSqr(pt - hitSurf.localPoints()[e.start()])
+ < matchTolerance
+ )
+ || (
+ magSqr(pt - hitSurf.localPoints()[e.end()])
+ < matchTolerance
+ )
+ )
+ {
+ continue;
+ }
+
+ nChanged++;
+
+ // Keep the points in the same place and move the edge
+ newPoints[hitSurfI].append(newPoints[surfI][pointI]);
+
+ // Move the points to the edges
+ //newPoints[pointI] = eHit.hitPoint();
+ //newPoints.append(eHit.hitPoint());
+
+ visitedFace[hitSurfI][faceI] = true;
+
+ // Split the other face.
+ greenRefine
+ (
+ hitSurf,
+ faceI,
+ eIndex,
+ newPoints[hitSurfI].size() - 1,
+ newFacesFromSplit
+ );
+
+ forAll(newFacesFromSplit, newFaceI)
+ {
+ if (newFaceI == 0)
+ {
+ newFaces[hitSurfI][faceI] = newFacesFromSplit[0];
+ }
+ else
+ {
+ newFaces[hitSurfI].append
+ (
+ newFacesFromSplit[newFaceI]
+ );
+ }
+ }
+ }
+ }
+ }
+
+ Info<< " Number of edges split = " << nChanged << endl;
+
+ forAll(newSurfaces, surfI)
+ {
+ newSurfaces.set
+ (
+ surfI,
+ new triSurfaceMesh
+ (
+ IOobject
+ (
+ "hookedSurface_" + surfs.names()[surfI],
+ runTime.constant(),
+ "triSurface",
+ runTime
+ ),
+ triSurface
+ (
+ newFaces[surfI],
+ newSurfaces[surfI].patches(),
+ pointField(newPoints[surfI])
+ )
+ )
+ );
+ }
+
+ } while (nChanged > 0);
+
+ Info<< endl;
+
+ forAll(newSurfaces, surfI)
+ {
+ const triSurfaceMesh& newSurf = newSurfaces[surfI];
+
+ Info<< "Writing hooked surface " << newSurf.searchableSurface::name()
+ << endl;
+
+ newSurf.searchableSurface::write();
+ }
+
+ Info<< "\nEnd\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceHookUp/surfaceHookUpDict b/applications/utilities/surface/surfaceHookUp/surfaceHookUpDict
new file mode 100644
index 0000000000000000000000000000000000000000..b25f21b5d8b64567b851650557d53af9beebaa96
--- /dev/null
+++ b/applications/utilities/surface/surfaceHookUp/surfaceHookUpDict
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object surfaceHookUpDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+surface1.stl {type triSurfaceMesh;}
+surface2.stl {type triSurfaceMesh;}
+
+
+// ************************************************************************* //
diff --git a/etc/config/settings.csh b/etc/config/settings.csh
index 2264c4ee03850cf37610384befd7b4238ba99bd2..a909e02de69799ddaa6c2a63f12d12b4a24fcb29 100644
--- a/etc/config/settings.csh
+++ b/etc/config/settings.csh
@@ -242,8 +242,7 @@ case ThirdParty:
# using clang - not gcc
setenv WM_CC 'clang'
setenv WM_CXX 'clang++'
- set clang_version=llvm-3.1
- #set clang_version=llvm-svn
+ set clang_version=llvm-3.2
breaksw
default:
echo
diff --git a/etc/config/settings.sh b/etc/config/settings.sh
index 5a77479b568f7e653a9fb4b0e81293f1357d6748..1d4800d471bd8462b3f4c89ca415f9142eeecd1d 100644
--- a/etc/config/settings.sh
+++ b/etc/config/settings.sh
@@ -262,8 +262,7 @@ OpenFOAM | ThirdParty)
# using clang - not gcc
export WM_CC='clang'
export WM_CXX='clang++'
- clang_version=llvm-3.1
- #clang_version=llvm-svn
+ clang_version=llvm-3.2
;;
*)
echo 1>&2
diff --git a/src/OpenFOAM/fields/Fields/Field/Field.C b/src/OpenFOAM/fields/Fields/Field/Field.C
index 8545c2aa9e3d85082938034eb8c6bbab82fc960e..334eba8e418962b2c39e6d4b114b09ce21736e7a 100644
--- a/src/OpenFOAM/fields/Fields/Field/Field.C
+++ b/src/OpenFOAM/fields/Fields/Field/Field.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) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -677,7 +677,6 @@ template<class Type>
template<class Form, class Cmpt, int nCmpt>
void Foam::Field<Type>::operator=(const VectorSpace<Form,Cmpt,nCmpt>& vs)
{
- typedef VectorSpace<Form,Cmpt,nCmpt> VSType;
TFOR_ALL_F_OP_S(Type, *this, =, VSType, vs)
}
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index ddf9564230faca1a7903b9f1929ae5d1dc3075b2..5f87a9ceef9f8e2de3b3df88a3b47455f57a5c5d 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -321,8 +321,6 @@ $(ddtSchemes)/CoEulerDdtScheme/CoEulerDdtSchemes.C
$(ddtSchemes)/SLTSDdtScheme/SLTSDdtSchemes.C
$(ddtSchemes)/localEulerDdtScheme/localEulerDdtSchemes.C
$(ddtSchemes)/backwardDdtScheme/backwardDdtSchemes.C
-$(ddtSchemes)/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C
-$(ddtSchemes)/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C
$(ddtSchemes)/CrankNicolsonDdtScheme/CrankNicolsonDdtSchemes.C
$(ddtSchemes)/boundedDdtScheme/boundedDdtSchemes.C
diff --git a/src/finiteVolume/cfdTools/general/include/alphaControls.H b/src/finiteVolume/cfdTools/general/include/alphaControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..b5ae1f94d83db4f45bd267e0a700a318f1482e79
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/include/alphaControls.H
@@ -0,0 +1,12 @@
+const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
+label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+Switch MULESCorr(alphaControls.lookup("MULESCorr"));
+
+if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
+{
+ FatalErrorIn(args.executable())
+ << "Sub-cycling alpha is only allowed for PISO operation, "
+ "i.e. when the number of outer-correctors = 1"
+ << exit(FatalError);
+}
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedPushedInternalValue/mappedFixedPushedInternalValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedPushedInternalValue/mappedFixedPushedInternalValueFvPatchField.C
index 2cea02cd9cd1d3fba125fff69cbf3227c9e6e69b..3550e414cf4b3203274f194e7325ed8f24e845c3 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedPushedInternalValue/mappedFixedPushedInternalValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedFixedPushedInternalValue/mappedFixedPushedInternalValueFvPatchField.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) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -95,8 +95,6 @@ mappedFixedPushedInternalValueFvPatchField
template<class Type>
void Foam::mappedFixedPushedInternalValueFvPatchField<Type>::updateCoeffs()
{
- typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
-
if (this->updated())
{
return;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C
deleted file mode 100644
index b46ca8e0f9959d36fd1661d1004340e7887583fa..0000000000000000000000000000000000000000
--- a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C
+++ /dev/null
@@ -1,708 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "boundedBackwardDdtScheme.H"
-#include "surfaceInterpolate.H"
-#include "fvcDiv.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-scalar boundedBackwardDdtScheme::deltaT_() const
-{
- return mesh().time().deltaTValue();
-}
-
-
-scalar boundedBackwardDdtScheme::deltaT0_() const
-{
- return mesh().time().deltaT0Value();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-tmp<volScalarField>
-boundedBackwardDdtScheme::fvcDdt
-(
- const dimensionedScalar& dt
-)
-{
- // No change compared to backward
-
- dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
-
- IOobject ddtIOobject
- (
- "ddt("+dt.name()+')',
- mesh().time().timeName(),
- mesh()
- );
-
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_();
-
- scalar coefft = 1 + deltaT/(deltaT + deltaT0);
- scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
- scalar coefft0 = coefft + coefft00;
-
- if (mesh().moving())
- {
- tmp<volScalarField> tdtdt
- (
- new volScalarField
- (
- ddtIOobject,
- mesh(),
- dimensionedScalar
- (
- "0",
- dt.dimensions()/dimTime,
- 0.0
- )
- )
- );
-
- tdtdt().internalField() = rDeltaT.value()*dt.value()*
- (
- coefft - (coefft0*mesh().V0() - coefft00*mesh().V00())/mesh().V()
- );
-
- return tdtdt;
- }
- else
- {
- return tmp<volScalarField>
- (
- new volScalarField
- (
- ddtIOobject,
- mesh(),
- dimensionedScalar
- (
- "0",
- dt.dimensions()/dimTime,
- 0.0
- ),
- calculatedFvPatchScalarField::typeName
- )
- );
- }
-}
-
-
-tmp<volScalarField>
-boundedBackwardDdtScheme::fvcDdt
-(
- const volScalarField& vf
-)
-{
- dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
-
- IOobject ddtIOobject
- (
- "ddt("+vf.name()+')',
- mesh().time().timeName(),
- mesh()
- );
-
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_(vf);
-
- // Calculate unboundedness indicator
- // Note: all times moved by one because access to internal field
- // copies current field into the old-time level.
- const volScalarField phict
- (
- mag
- (
- vf.oldTime().oldTime()
- - vf.oldTime().oldTime().oldTime()
- )/
- (
- mag
- (
- vf.oldTime()
- - vf.oldTime().oldTime()
- )
- + dimensionedScalar("small", vf.dimensions(), VSMALL)
- )
- );
-
- const volScalarField limiter(pos(phict) - pos(phict - scalar(1)));
- const volScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
- const volScalarField coefft00
- (
- limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
- );
- const volScalarField coefft0(coefft + coefft00);
-
- if (mesh().moving())
- {
- return tmp<volScalarField>
- (
- new volScalarField
- (
- ddtIOobject,
- mesh(),
- rDeltaT.dimensions()*vf.dimensions(),
- rDeltaT.value()*
- (
- coefft*vf.internalField() -
- (
- coefft0.internalField()
- *vf.oldTime().internalField()*mesh().V0()
- - coefft00.internalField()
- *vf.oldTime().oldTime().internalField()
- *mesh().V00()
- )/mesh().V()
- ),
- rDeltaT.value()*
- (
- coefft.boundaryField()*vf.boundaryField() -
- (
- coefft0.boundaryField()*
- vf.oldTime().boundaryField()
- - coefft00.boundaryField()*
- vf.oldTime().oldTime().boundaryField()
- )
- )
- )
- );
- }
- else
- {
- return tmp<volScalarField>
- (
- new volScalarField
- (
- ddtIOobject,
- rDeltaT*
- (
- coefft*vf
- - coefft0*vf.oldTime()
- + coefft00*vf.oldTime().oldTime()
- )
- )
- );
- }
-}
-
-
-tmp<volScalarField>
-boundedBackwardDdtScheme::fvcDdt
-(
- const dimensionedScalar& rho,
- const volScalarField& vf
-)
-{
- dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
-
- IOobject ddtIOobject
- (
- "ddt("+rho.name()+','+vf.name()+')',
- mesh().time().timeName(),
- mesh()
- );
-
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_(vf);
-
- // Calculate unboundedness indicator
- // Note: all times moved by one because access to internal field
- // copies current field into the old-time level.
- const volScalarField phict
- (
- mag
- (
- vf.oldTime().oldTime()
- - vf.oldTime().oldTime().oldTime()
- )/
- (
- mag
- (
- vf.oldTime()
- - vf.oldTime().oldTime()
- )
- + dimensionedScalar("small", vf.dimensions(), VSMALL)
- )
- );
-
- const volScalarField limiter(pos(phict) - pos(phict - scalar(1)));
- const volScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
- const volScalarField coefft00
- (
- limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
- );
- const volScalarField coefft0(coefft + coefft00);
-
- if (mesh().moving())
- {
- return tmp<volScalarField>
- (
- new volScalarField
- (
- ddtIOobject,
- mesh(),
- rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
- rDeltaT.value()*rho.value()*
- (
- coefft*vf.internalField() -
- (
- coefft0.internalField()*
- vf.oldTime().internalField()*mesh().V0()
- - coefft00.internalField()*
- vf.oldTime().oldTime().internalField()
- *mesh().V00()
- )/mesh().V()
- ),
- rDeltaT.value()*rho.value()*
- (
- coefft.boundaryField()*vf.boundaryField() -
- (
- coefft0.boundaryField()*
- vf.oldTime().boundaryField()
- - coefft00.boundaryField()*
- vf.oldTime().oldTime().boundaryField()
- )
- )
- )
- );
- }
- else
- {
- return tmp<volScalarField>
- (
- new volScalarField
- (
- ddtIOobject,
- rDeltaT*rho*
- (
- coefft*vf
- - coefft0*vf.oldTime()
- + coefft00*vf.oldTime().oldTime()
- )
- )
- );
- }
-}
-
-
-tmp<volScalarField>
-boundedBackwardDdtScheme::fvcDdt
-(
- const volScalarField& rho,
- const volScalarField& vf
-)
-{
- dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
-
- IOobject ddtIOobject
- (
- "ddt("+rho.name()+','+vf.name()+')',
- mesh().time().timeName(),
- mesh()
- );
-
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_(vf);
-
- // Calculate unboundedness indicator
- // Note: all times moved by one because access to internal field
- // copies current field into the old-time level.
- const volScalarField phict
- (
- mag
- (
- rho.oldTime().oldTime()*vf.oldTime().oldTime()
- - rho.oldTime().oldTime().oldTime()*vf.oldTime().oldTime().oldTime()
- )/
- (
- mag
- (
- rho.oldTime()*vf.oldTime()
- - rho.oldTime().oldTime()*vf.oldTime().oldTime()
- )
- + dimensionedScalar("small", rho.dimensions()*vf.dimensions(), VSMALL)
- )
- );
-
- const volScalarField limiter(pos(phict) - pos(phict - scalar(1)));
- const volScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
- const volScalarField coefft00
- (
- limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
- );
- const volScalarField coefft0(coefft + coefft00);
-
- if (mesh().moving())
- {
- return tmp<volScalarField>
- (
- new volScalarField
- (
- ddtIOobject,
- mesh(),
- rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
- rDeltaT.value()*
- (
- coefft*rho.internalField()*vf.internalField() -
- (
- coefft0.internalField()*
- rho.oldTime().internalField()*
- vf.oldTime().internalField()*mesh().V0()
- - coefft00.internalField()*
- rho.oldTime().oldTime().internalField()
- *vf.oldTime().oldTime().internalField()*mesh().V00()
- )/mesh().V()
- ),
- rDeltaT.value()*
- (
- coefft.boundaryField()*vf.boundaryField() -
- (
- coefft0.boundaryField()*
- rho.oldTime().boundaryField()*
- vf.oldTime().boundaryField()
- - coefft00.boundaryField()*
- rho.oldTime().oldTime().boundaryField()*
- vf.oldTime().oldTime().boundaryField()
- )
- )
- )
- );
- }
- else
- {
- return tmp<volScalarField>
- (
- new volScalarField
- (
- ddtIOobject,
- rDeltaT*
- (
- coefft*rho*vf
- - coefft0*rho.oldTime()*vf.oldTime()
- + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
- )
- )
- );
- }
-}
-
-
-tmp<fvScalarMatrix>
-boundedBackwardDdtScheme::fvmDdt
-(
- const volScalarField& vf
-)
-{
- tmp<fvScalarMatrix> tfvm
- (
- new fvScalarMatrix
- (
- vf,
- vf.dimensions()*dimVol/dimTime
- )
- );
-
- fvScalarMatrix& fvm = tfvm();
-
- scalar rDeltaT = 1.0/deltaT_();
-
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_(vf);
-
- // Calculate unboundedness indicator
- // Note: all times moved by one because access to internal field
- // copies current field into the old-time level.
- const scalarField phict
- (
- mag
- (
- vf.oldTime().oldTime().internalField()
- - vf.oldTime().oldTime().oldTime().internalField()
- )/
- (
- mag
- (
- vf.oldTime().internalField()
- - vf.oldTime().oldTime().internalField()
- )
- + VSMALL
- )
- );
-
- const scalarField limiter(pos(phict) - pos(phict - 1.0));
- const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
- const scalarField coefft00
- (
- limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
- );
- const scalarField coefft0(coefft + coefft00);
-
- fvm.diag() = (coefft*rDeltaT)*mesh().V();
-
- if (mesh().moving())
- {
- fvm.source() = rDeltaT*
- (
- coefft0*vf.oldTime().internalField()*mesh().V0()
- - coefft00*vf.oldTime().oldTime().internalField()
- *mesh().V00()
- );
- }
- else
- {
- fvm.source() = rDeltaT*mesh().V()*
- (
- coefft0*vf.oldTime().internalField()
- - coefft00*vf.oldTime().oldTime().internalField()
- );
- }
-
- return tfvm;
-}
-
-
-tmp<fvScalarMatrix>
-boundedBackwardDdtScheme::fvmDdt
-(
- const dimensionedScalar& rho,
- const volScalarField& vf
-)
-{
- tmp<fvScalarMatrix> tfvm
- (
- new fvScalarMatrix
- (
- vf,
- rho.dimensions()*vf.dimensions()*dimVol/dimTime
- )
- );
- fvScalarMatrix& fvm = tfvm();
-
- scalar rDeltaT = 1.0/deltaT_();
-
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_(vf);
-
- // Calculate unboundedness indicator
- // Note: all times moved by one because access to internal field
- // copies current field into the old-time level.
- const scalarField phict
- (
- mag
- (
- vf.oldTime().oldTime().internalField()
- - vf.oldTime().oldTime().oldTime().internalField()
- )/
- (
- mag
- (
- vf.oldTime().internalField()
- - vf.oldTime().oldTime().internalField()
- )
- + VSMALL
- )
- );
-
- const scalarField limiter(pos(phict) - pos(phict - 1.0));
- const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
- const scalarField coefft00
- (
- limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
- );
- const scalarField coefft0(coefft + coefft00);
-
- fvm.diag() = (coefft*rDeltaT*rho.value())*mesh().V();
-
- if (mesh().moving())
- {
- fvm.source() = rDeltaT*rho.value()*
- (
- coefft0*vf.oldTime().internalField()*mesh().V0()
- - coefft00*vf.oldTime().oldTime().internalField()
- *mesh().V00()
- );
- }
- else
- {
- fvm.source() = rDeltaT*mesh().V()*rho.value()*
- (
- coefft0*vf.oldTime().internalField()
- - coefft00*vf.oldTime().oldTime().internalField()
- );
- }
-
- return tfvm;
-}
-
-
-tmp<fvScalarMatrix>
-boundedBackwardDdtScheme::fvmDdt
-(
- const volScalarField& rho,
- const volScalarField& vf
-)
-{
- tmp<fvScalarMatrix> tfvm
- (
- new fvScalarMatrix
- (
- vf,
- rho.dimensions()*vf.dimensions()*dimVol/dimTime
- )
- );
- fvScalarMatrix& fvm = tfvm();
-
- scalar rDeltaT = 1.0/deltaT_();
-
- scalar deltaT = deltaT_();
- scalar deltaT0 = deltaT0_(vf);
-
- // Calculate unboundedness indicator
- // Note: all times moved by one because access to internal field
- // copies current field into the old-time level.
- const scalarField phict
- (
- mag
- (
- rho.oldTime().oldTime().internalField()*
- vf.oldTime().oldTime().internalField()
- - rho.oldTime().oldTime().oldTime().internalField()*
- vf.oldTime().oldTime().oldTime().internalField()
- )/
- (
- mag
- (
- rho.oldTime().internalField()*
- vf.oldTime().internalField()
- - rho.oldTime().oldTime().internalField()*
- vf.oldTime().oldTime().internalField()
- )
- + VSMALL
- )
- );
-
- const scalarField limiter(pos(phict) - pos(phict - 1.0));
- const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
- const scalarField coefft00
- (
- limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
- );
- const scalarField coefft0(coefft + coefft00);
-
- fvm.diag() = (coefft*rDeltaT)*rho.internalField()*mesh().V();
-
- if (mesh().moving())
- {
- fvm.source() = rDeltaT*
- (
- coefft0*rho.oldTime().internalField()
- *vf.oldTime().internalField()*mesh().V0()
- - coefft00*rho.oldTime().oldTime().internalField()
- *vf.oldTime().oldTime().internalField()*mesh().V00()
- );
- }
- else
- {
- fvm.source() = rDeltaT*mesh().V()*
- (
- coefft0*rho.oldTime().internalField()
- *vf.oldTime().internalField()
- - coefft00*rho.oldTime().oldTime().internalField()
- *vf.oldTime().oldTime().internalField()
- );
- }
-
- return tfvm;
-}
-
-
-tmp<surfaceScalarField> boundedBackwardDdtScheme::fvcDdtPhiCorr
-(
- const volScalarField& rA,
- const volScalarField& U,
- const surfaceScalarField& phi
-)
-{
- notImplemented
- (
- "boundedBackwardDdtScheme::fvcDdtPhiCorr"
- );
-
- return surfaceScalarField::null();
-}
-
-
-tmp<surfaceScalarField> boundedBackwardDdtScheme::fvcDdtPhiCorr
-(
- const volScalarField& rA,
- const volScalarField& rho,
- const volScalarField& U,
- const surfaceScalarField& phi
-)
-{
- notImplemented
- (
- "boundedBackwardDdtScheme::fvcDdtPhiCorr"
- );
-
- return surfaceScalarField::null();
-}
-
-
-tmp<surfaceScalarField> boundedBackwardDdtScheme::meshPhi
-(
- const volScalarField& vf
-)
-{
- notImplemented
- (
- "boundedBackwardDdtScheme::meshPhi(const volScalarField& vf)"
- );
-
- return surfaceScalarField::null();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.H
deleted file mode 100644
index 900e2a7993295db1780dd4caf107d7699e69eaf2..0000000000000000000000000000000000000000
--- a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.H
+++ /dev/null
@@ -1,193 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
-
-Class
- Foam::fv::boundedBackwardDdtScheme
-
-Description
- Second-order bounded-backward-differencing ddt using the current and
- two previous time-step values.
-
-SourceFiles
- boundedBackwardDdtScheme.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef boundedBackwardDdtScheme_H
-#define boundedBackwardDdtScheme_H
-
-#include "ddtScheme.H"
-#include "fvMatrices.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-/*---------------------------------------------------------------------------*\
- Class boundedBackwardDdtScheme Declaration
-\*---------------------------------------------------------------------------*/
-
-class boundedBackwardDdtScheme
-:
- public fv::ddtScheme<scalar>
-{
- // Private Member Functions
-
- //- Return the current time-step
- scalar deltaT_() const;
-
- //- Return the previous time-step
- scalar deltaT0_() const;
-
- //- Return the previous time-step or GREAT if the old timestep field
- // wasn't available in which case Euler ddt is used
- template<class GeoField>
- scalar deltaT0_(const GeoField& vf) const
- {
- if (vf.oldTime().timeIndex() == vf.oldTime().oldTime().timeIndex())
- {
- return GREAT;
- }
- else
- {
- return deltaT0_();
- }
- }
-
-
- //- Disallow default bitwise copy construct
- boundedBackwardDdtScheme(const boundedBackwardDdtScheme&);
-
- //- Disallow default bitwise assignment
- void operator=(const boundedBackwardDdtScheme&);
-
-
-public:
-
- //- Runtime type information
- TypeName("boundedBackward");
-
-
- // Constructors
-
- //- Construct from mesh
- boundedBackwardDdtScheme(const fvMesh& mesh)
- :
- ddtScheme<scalar>(mesh)
- {}
-
- //- Construct from mesh and Istream
- boundedBackwardDdtScheme(const fvMesh& mesh, Istream& is)
- :
- ddtScheme<scalar>(mesh, is)
- {}
-
-
- // Member Functions
-
- //- Return mesh reference
- const fvMesh& mesh() const
- {
- return fv::ddtScheme<scalar>::mesh();
- }
-
- tmp<volScalarField> fvcDdt
- (
- const dimensionedScalar&
- );
-
- tmp<volScalarField> fvcDdt
- (
- const volScalarField&
- );
-
- tmp<volScalarField> fvcDdt
- (
- const dimensionedScalar&,
- const volScalarField&
- );
-
- tmp<volScalarField> fvcDdt
- (
- const volScalarField&,
- const volScalarField&
- );
-
- tmp<fvScalarMatrix> fvmDdt
- (
- const volScalarField&
- );
-
- tmp<fvScalarMatrix> fvmDdt
- (
- const dimensionedScalar&,
- const volScalarField&
- );
-
- tmp<fvScalarMatrix> fvmDdt
- (
- const volScalarField&,
- const volScalarField&
- );
-
- tmp<surfaceScalarField> fvcDdtPhiCorr
- (
- const volScalarField& rA,
- const volScalarField& U,
- const surfaceScalarField& phi
- );
-
- tmp<surfaceScalarField> fvcDdtPhiCorr
- (
- const volScalarField& rA,
- const volScalarField& rho,
- const volScalarField& U,
- const surfaceScalarField& phi
- );
-
- tmp<surfaceScalarField> meshPhi
- (
- const volScalarField&
- );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C
deleted file mode 100644
index 1504b8ba3b39dc7b74f5c070ceb99c20824cbc88..0000000000000000000000000000000000000000
--- a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C
+++ /dev/null
@@ -1,44 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "boundedBackwardDdtScheme.H"
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace fv
-{
-
-defineTypeNameAndDebug(boundedBackwardDdtScheme, 0);
-
-ddtScheme<scalar>::addIstreamConstructorToTable<boundedBackwardDdtScheme>
- addboundedBackwardDdtSchemeIstreamConstructorToTable_;
-
-}
-}
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
index 62ca855e785d375bf752f19119b18d667c554117..ff5994e670ce942784ddbdb283ebc4247896a603 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
@@ -48,6 +48,25 @@ void Foam::MULES::explicitSolve
}
+void Foam::MULES::correct
+(
+ volScalarField& psi,
+ surfaceScalarField& phiPsiCorr,
+ const scalar psiMax,
+ const scalar psiMin
+)
+{
+ correct
+ (
+ geometricOneField(),
+ psi,
+ phiPsiCorr,
+ zeroField(), zeroField(),
+ psiMax, psiMin
+ );
+}
+
+
void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
{
forAll(phiPsiCorrs[0], facei)
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
index 2307c289e4007bddb4436326d04db4049d453d32..c1d824b9503a16fce8f27befc168e3a40b754e66 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
@@ -92,6 +92,38 @@ 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
+(
+ const RhoType& rho,
+ volScalarField& psi,
+ surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin
+);
+
+void correct
+(
+ volScalarField& psi,
+ surfaceScalarField& phiCorr,
+ const scalar psiMax,
+ const scalar psiMin
+);
+
+
template<class RhoType, class SpType, class SuType>
void limiter
(
@@ -122,6 +154,35 @@ void limit
const bool returnCorr
);
+
+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>
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index a3cb9f17e9d40e52a7b5c850975eb00dbcb27eed..e0718f923a91b2d3f26fd7db813447c8f9822f8d 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -96,6 +96,74 @@ void Foam::MULES::explicitSolve
}
+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
+)
+{
+ 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());
+ }
+
+ psi.correctBoundaryConditions();
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::correct
+(
+ const RhoType& rho,
+ volScalarField& psi,
+ surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin
+)
+{
+ const fvMesh& mesh = psi.mesh();
+
+ const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+ label nLimiterIter
+ (
+ readLabel(MULEScontrols.lookup("nLimiterIter"))
+ );
+
+ limitCorr(rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+ correct(rho, psi, phiCorr, Sp, Su);
+}
+
+
template<class RhoType, class SpType, class SuType>
void Foam::MULES::limiter
(
@@ -511,6 +579,364 @@ 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] > 0)
+ {
+ 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/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
index 48d73efb2e316a6d00bbc44893a7a6aa0206c460..9a7462127cded0319c428f707a19272ccf94e150 100644
--- a/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.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) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -43,8 +43,10 @@ bool Foam::pairPatchAgglomeration::continueAgglomerating
)
{
// Check the need for further agglomeration on all processors
- bool contAgg = nCoarseFaces >= nFacesInCoarsestLevel_;
- reduce(contAgg, andOp<bool>());
+ label localnCoarseFaces = nCoarseFaces;
+// reduce(localnCoarseFaces, sumOp<label>());
+ bool contAgg = localnCoarseFaces >= nFacesInCoarsestLevel_;
+ //reduce(contAgg, andOp<bool>());
return contAgg;
}
@@ -263,6 +265,11 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
) << "min(fineToCoarse) == -1" << exit(FatalError);
}
+ if (fineToCoarse.size() == 0)
+ {
+ return true;
+ }
+
if (fineToCoarse.size() != patch.size())
{
FatalErrorIn
@@ -282,6 +289,7 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
const label nCoarseI = max(fineToCoarse)+1;
List<face> patchFaces(nCoarseI);
+
// Patch faces per agglomeration
labelListList coarseToFine(invertOneToMany(nCoarseI, fineToCoarse));
@@ -335,6 +343,7 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
patch.points()
)
);
+
return true;
}
@@ -343,59 +352,73 @@ void Foam::pairPatchAgglomeration:: agglomerate()
{
label nPairLevels = 0;
label nCreatedLevels = 1; //0 level is the base patch
+ label nCoarseFaces = 0;
+ label nCoarseFacesOld = 0;
while (nCreatedLevels < maxLevels_)
{
- label nCoarseCells = -1;
-
const bPatch& patch = patchLevels_[nCreatedLevels - 1];
tmp<labelField> finalAgglomPtr(new labelField(patch.size()));
bool agglomOK = false;
- while (!agglomOK && patch.size())
+ while (!agglomOK)
{
finalAgglomPtr = agglomerateOneLevel
(
- nCoarseCells,
+ nCoarseFaces,
patch
);
- agglomOK = agglomeratePatch
- (
- patch,
- finalAgglomPtr,
- nCreatedLevels
- );
- }
+ if (nCoarseFaces > 0)
+ {
+ agglomOK = agglomeratePatch
+ (
+ patch,
+ finalAgglomPtr,
+ nCreatedLevels
+ );
- nFaces_[nCreatedLevels] = nCoarseCells;
- restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
- mapBaseToTopAgglom(nCreatedLevels);
+ restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
+ mapBaseToTopAgglom(nCreatedLevels);
+ setEdgeWeights(nCreatedLevels);
- if (!continueAgglomerating(nCoarseCells))
- {
- break;
+ if (nPairLevels % mergeLevels_)
+ {
+ combineLevels(nCreatedLevels);
+ }
+ else
+ {
+ nCreatedLevels++;
+ }
+
+ nPairLevels++;
+ }
+ else
+ {
+ agglomOK = true;
+ }
+ reduce(nCoarseFaces, sumOp<label>());
}
- setEdgeWeights(nCreatedLevels);
+ nFaces_[nCreatedLevels] = nCoarseFaces;
- if (nPairLevels % mergeLevels_)
- {
- combineLevels(nCreatedLevels);
- }
- else
+ if
+ (
+ !continueAgglomerating(nCoarseFaces)
+ || (nCoarseFacesOld == nCoarseFaces)
+ )
{
- nCreatedLevels++;
+ break;
}
- nPairLevels++;
+ nCoarseFacesOld = nCoarseFaces;
}
}
Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
(
- label& nCoarseCells,
+ label& nCoarseFaces,
const bPatch& patch
)
{
@@ -406,7 +429,7 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
const labelListList& faceFaces = patch.faceFaces();
- nCoarseCells = 0;
+ nCoarseFaces = 0;
forAll(faceFaces, facei)
{
@@ -440,9 +463,9 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
if (matchFaceNo >= 0)
{
// Make a new group
- coarseCellMap[matchFaceNo] = nCoarseCells;
- coarseCellMap[matchFaceNeibNo] = nCoarseCells;
- nCoarseCells++;
+ coarseCellMap[matchFaceNo] = nCoarseFaces;
+ coarseCellMap[matchFaceNeibNo] = nCoarseFaces;
+ nCoarseFaces++;
}
else
{
@@ -475,8 +498,8 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
else
{
// if not create single-cell "clusters" for each
- coarseCellMap[facei] = nCoarseCells;
- nCoarseCells ++;
+ coarseCellMap[facei] = nCoarseFaces;
+ nCoarseFaces ++;
}
}
}
diff --git a/src/meshTools/indexedOctree/treeDataEdge.C b/src/meshTools/indexedOctree/treeDataEdge.C
index 4d0b9b4d813624bb635424c125b1f61e7e3f0d07..ffd0969e8f5612ab5d838bf0a7b851c137e1761b 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.C
+++ b/src/meshTools/indexedOctree/treeDataEdge.C
@@ -105,6 +105,17 @@ Foam::treeDataEdge::findNearestOp::findNearestOp
{}
+Foam::treeDataEdge::findNearestOpSubset::findNearestOpSubset
+(
+ const indexedOctree<treeDataEdge>& tree,
+ DynamicList<label>& shapeMask
+)
+:
+ tree_(tree),
+ shapeMask_(shapeMask)
+{}
+
+
Foam::treeDataEdge::findIntersectOp::findIntersectOp
(
const indexedOctree<treeDataEdge>& tree
@@ -270,6 +281,105 @@ void Foam::treeDataEdge::findNearestOp::operator()
}
+void Foam::treeDataEdge::findNearestOpSubset::operator()
+(
+ const labelUList& indices,
+ const point& sample,
+
+ scalar& nearestDistSqr,
+ label& minIndex,
+ point& nearestPoint
+) const
+{
+ const treeDataEdge& shape = tree_.shapes();
+
+ forAll(indices, i)
+ {
+ const label index = indices[i];
+ const label edgeIndex = shape.edgeLabels()[index];
+
+ if (!shapeMask_.empty() && findIndex(shapeMask_, edgeIndex) != -1)
+ {
+ continue;
+ }
+
+ const edge& e = shape.edges()[edgeIndex];
+
+ pointHit nearHit = e.line(shape.points()).nearestDist(sample);
+
+ scalar distSqr = sqr(nearHit.distance());
+
+ if (distSqr < nearestDistSqr)
+ {
+ nearestDistSqr = distSqr;
+ minIndex = index;
+ nearestPoint = nearHit.rawPoint();
+ }
+ }
+}
+
+
+void Foam::treeDataEdge::findNearestOpSubset::operator()
+(
+ const labelUList& indices,
+ const linePointRef& ln,
+
+ treeBoundBox& tightest,
+ label& minIndex,
+ point& linePoint,
+ point& nearestPoint
+) const
+{
+ const treeDataEdge& shape = tree_.shapes();
+
+ // Best so far
+ scalar nearestDistSqr = magSqr(linePoint - nearestPoint);
+
+ forAll(indices, i)
+ {
+ const label index = indices[i];
+ const label edgeIndex = shape.edgeLabels()[index];
+
+ if (!shapeMask_.empty() && findIndex(shapeMask_, edgeIndex) != -1)
+ {
+ continue;
+ }
+
+ const edge& e = shape.edges()[edgeIndex];
+
+ // Note: could do bb test ? Worthwhile?
+
+ // Nearest point on line
+ point ePoint, lnPt;
+ scalar dist = e.line(shape.points()).nearestDist(ln, ePoint, lnPt);
+ scalar distSqr = sqr(dist);
+
+ if (distSqr < nearestDistSqr)
+ {
+ nearestDistSqr = distSqr;
+ minIndex = index;
+ linePoint = lnPt;
+ nearestPoint = ePoint;
+
+ {
+ point& minPt = tightest.min();
+ minPt = min(ln.start(), ln.end());
+ minPt.x() -= dist;
+ minPt.y() -= dist;
+ minPt.z() -= dist;
+ }
+ {
+ point& maxPt = tightest.max();
+ maxPt = max(ln.start(), ln.end());
+ maxPt.x() += dist;
+ maxPt.y() += dist;
+ maxPt.z() += dist;
+ }
+ }
+ }
+}
+
+
bool Foam::treeDataEdge::findIntersectOp::operator()
(
const label index,
diff --git a/src/meshTools/indexedOctree/treeDataEdge.H b/src/meshTools/indexedOctree/treeDataEdge.H
index 12fd0d6aa10ba98c5db55977549b4775d58ca4d7..5a804d5f8e337228785fa3dc55920d9400409f2c 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.H
+++ b/src/meshTools/indexedOctree/treeDataEdge.H
@@ -118,6 +118,43 @@ public:
) const;
};
+ class findNearestOpSubset
+ {
+ const indexedOctree<treeDataEdge>& tree_;
+
+ DynamicList<label>& shapeMask_;
+
+ public:
+
+ findNearestOpSubset
+ (
+ const indexedOctree<treeDataEdge>& tree,
+ DynamicList<label>& shapeMask
+ );
+
+ void operator()
+ (
+ const labelUList& indices,
+ const point& sample,
+
+ scalar& nearestDistSqr,
+ label& minIndex,
+ point& nearestPoint
+ ) const;
+
+ //- Calculates nearest (to line) point in shape.
+ // Returns point and distance (squared)
+ void operator()
+ (
+ const labelUList& indices,
+ const linePointRef& ln,
+
+ treeBoundBox& tightest,
+ label& minIndex,
+ point& linePoint,
+ point& nearestPoint
+ ) const;
+ };
class findIntersectOp
{
diff --git a/src/postProcessing/foamCalcFunctions/basic/addSubtract/addSubtract.C b/src/postProcessing/foamCalcFunctions/basic/addSubtract/addSubtract.C
index c284e2c0426ad64411726c2f67e8a6dcf7be2240..a4a7941055f00c0979cdac95bdfc29cbafa73a06 100644
--- a/src/postProcessing/foamCalcFunctions/basic/addSubtract/addSubtract.C
+++ b/src/postProcessing/foamCalcFunctions/basic/addSubtract/addSubtract.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) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -297,4 +297,3 @@ void Foam::calcTypes::addSubtract::calc
// ************************************************************************* //
-
diff --git a/src/postProcessing/foamCalcFunctions/basic/addSubtract/writeAddSubtractField.C b/src/postProcessing/foamCalcFunctions/basic/addSubtract/writeAddSubtractField.C
index e069cd1fc08abf8856c3cc038d76eaad704c816d..8304569ae4a65ed496f242a024158e7eccafe244 100644
--- a/src/postProcessing/foamCalcFunctions/basic/addSubtract/writeAddSubtractField.C
+++ b/src/postProcessing/foamCalcFunctions/basic/addSubtract/writeAddSubtractField.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) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -72,7 +72,9 @@ void Foam::calcTypes::addSubtract::writeAddSubtractField
mesh,
IOobject::NO_READ
),
- calcMode_ == ADD ? baseField + addField : baseField - addField
+ calcMode_ == ADD
+ ? (baseField + addField)()
+ : (baseField - addField)()
);
newField.write();
}
diff --git a/src/transportModels/interfaceProperties/interfaceProperties.C b/src/transportModels/interfaceProperties/interfaceProperties.C
index 63a1340ee92ddf3fa3dab052365f2e35bcfa3ee4..cfa20630aed40bb294fb050a4e4091bb61b41b39 100644
--- a/src/transportModels/interfaceProperties/interfaceProperties.C
+++ b/src/transportModels/interfaceProperties/interfaceProperties.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) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -159,7 +159,7 @@ Foam::interfaceProperties::interfaceProperties
(
readScalar
(
- alpha1.mesh().solutionDict().subDict("PIMPLE").lookup("cAlpha")
+ alpha1.mesh().solverDict(alpha1.name()).lookup("cAlpha")
)
),
sigma_(dict.lookup("sigma")),
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict
index f190eaf2ff82595058739bfa82a38aed4bc8ba1f..d1a8467089474642d52441dad9ba86ce154c204c 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict
@@ -14,37 +14,37 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-application interPhaseChangeFoam;
+application interPhaseChangeFoam;
-startFrom latestTime;
+startFrom latestTime;
-startTime 0;
+startTime 0;
-stopAt endTime;
+stopAt endTime;
-endTime 0.05;
+endTime 0.05;
-deltaT 1e-8;
+deltaT 1e-8;
-writeControl adjustableRunTime;
+writeControl adjustableRunTime;
-writeInterval 0.001;
+writeInterval 0.001;
-purgeWrite 0;
+purgeWrite 0;
-writeFormat ascii;
+writeFormat ascii;
-writePrecision 6;
+writePrecision 6;
-writeCompression uncompressed;
+writeCompression uncompressed;
-timeFormat general;
+timeFormat general;
-runTimeModifiable yes;
+runTimeModifiable yes;
-adjustTimeStep on;
+adjustTimeStep on;
-maxCo 1.0;
+maxCo 5;
// ************************************************************************* //
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
index 9143f603f4b73ec4b100fd003d242e65af899c57..8ddd1fe903d3c4bd9293c1afe4e9dbafb8a8ffd5 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
@@ -32,6 +32,7 @@ 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;
}
@@ -47,10 +48,7 @@ laplacianSchemes
snGradSchemes
{
- default none;
- snGrad(pd) limited corrected 0.5;
- snGrad(rho) limited corrected 0.5;
- snGrad(alpha1) limited corrected 0.5;
+ default limited corrected 0.5;
}
fluxRequired
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution
index bfbd5df2715b0f33a0bb2b07e93fcc83523c550a..a1f8f8ffd5e3c2376a31e3a25f9f49fe036cfd3d 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution
@@ -18,22 +18,24 @@ solvers
{
alpha1
{
- maxUnboundedness 1e-5;
- CoCoeff 2;
- maxIter 5;
- nLimiterIter 2;
+ cAlpha 0;
+ nAlphaCorr 2;
+ nAlphaSubCycles 1;
- solver PBiCG;
- preconditioner DILU;
- tolerance 1e-12;
- relTol 0.1;
+ MULESCorr yes;
+ nLimiterIter 5;
+
+ solver PBiCG;
+ preconditioner DILU;
+ tolerance 1e-8;
+ relTol 0;
};
- U
+ "U.*"
{
solver PBiCG;
preconditioner DILU;
- tolerance 1e-06;
+ tolerance 1e-6;
relTol 0;
};
@@ -96,10 +98,6 @@ PIMPLE
nOuterCorrectors 1;
nCorrectors 3;
nNonOrthogonalCorrectors 0;
-
- cAlpha 0;
- nAlphaCorr 1;
- nAlphaSubCycles 1;
}