How to use PIMPLE properly?

floquation · February 14, 2017, 09:16

Background:
I am simulating a rising bubble in OF41 with the interFoam solver.
To do so, I am currently using a uniform mesh with cubic cells (aspect ratio 1).
There are a total of 90x90x180 cells, which corresponds to 15 cells / bubble diameter. (This is relatively coarse, but should be fine enough to test the below.)

This case converges pretty well using only one outer (SIMPLE) iteration and just three inner (PISO) iterations with maxCo=0.5.
I decided to check whether using the PIMPLE algorithm - following Tobi's guide - might speed up my simulation.

Problem (I mean, challenge):
The PIMPLE algorithm slows down my simulation considerably, which is because the solution does not (barely) improve with successive outer iterations.
maxCo is only 2. (Which means that each timestep must at most take a factor four longer to be faster than the maxCo=0.5 simulation.)

Log file (snippet):

Code:

Courant Number mean: 0.000976416 max: 1.44982
Interface Courant Number mean: 0.000135741 max: 1.05712
deltaT = 7.14286e-05
Time = 0.0035

PIMPLE: iteration 1
smoothSolver:  Solving for alpha.water, Initial residual = 0.00797562, Final residual = 6.12225e-08, No Iterations 3
Phase-1 volume fraction = 0.998807  Min(alpha.water) = 2.79909e-07  Max(alpha.water) = 1.0002
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.998807  Min(alpha.water) = -0.000994  Max(alpha.water) = 1.0002
GAMG:  Solving for p_rgh, Initial residual = 0.135337, Final residual = 0.000314987, No Iterations 2
time step continuity errors : sum local = 5.46316e-07, global = 1.47295e-10, cumulative = -2.1714e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0940134, Final residual = 0.000266503, No Iterations 2
time step continuity errors : sum local = 4.6467e-07, global = -6.43504e-10, cumulative = -2.17783e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0941178, Final residual = 0.00025847, No Iterations 2
time step continuity errors : sum local = 4.50745e-07, global = -8.69603e-10, cumulative = -2.18653e-07
PIMPLE: iteration 2
smoothSolver:  Solving for alpha.water, Initial residual = 0.00645469, Final residual = 1.00334e-09, No Iterations 3
Phase-1 volume fraction = 0.998807  Min(alpha.water) = 2.06368e-07  Max(alpha.water) = 1.00028
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.998807  Min(alpha.water) = -0.00122852  Max(alpha.water) = 1.00028
GAMG:  Solving for p_rgh, Initial residual = 0.086505, Final residual = 0.000338267, No Iterations 2
time step continuity errors : sum local = 7.40111e-07, global = -3.00883e-10, cumulative = -2.18954e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0606414, Final residual = 0.000265495, No Iterations 2
time step continuity errors : sum local = 5.8013e-07, global = -7.78486e-10, cumulative = -2.19732e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0606601, Final residual = 0.000258879, No Iterations 2
time step continuity errors : sum local = 5.65676e-07, global = -8.99989e-10, cumulative = -2.20632e-07
PIMPLE: iteration 3
smoothSolver:  Solving for alpha.water, Initial residual = 0.00582153, Final residual = 2.94448e-09, No Iterations 3
Phase-1 volume fraction = 0.998807  Min(alpha.water) = 2.78524e-07  Max(alpha.water) = 1.00023
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.998807  Min(alpha.water) = -0.00139389  Max(alpha.water) = 1.00023
GAMG:  Solving for p_rgh, Initial residual = 0.080493, Final residual = 0.000340901, No Iterations 2
time step continuity errors : sum local = 6.77737e-07, global = -1.80287e-10, cumulative = -2.20813e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0560016, Final residual = 0.000246158, No Iterations 2
time step continuity errors : sum local = 4.90191e-07, global = -1.04078e-09, cumulative = -2.21853e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0560882, Final residual = 0.000241618, No Iterations 2
time step continuity errors : sum local = 4.81123e-07, global = -1.22133e-09, cumulative = -2.23075e-07
PIMPLE: iteration 4
smoothSolver:  Solving for alpha.water, Initial residual = 0.00559932, Final residual = 9.83465e-10, No Iterations 3
Phase-1 volume fraction = 0.998807  Min(alpha.water) = 2.60849e-07  Max(alpha.water) = 1.0002
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.998807  Min(alpha.water) = -0.000778589  Max(alpha.water) = 1.0002
GAMG:  Solving for p_rgh, Initial residual = 0.0612285, Final residual = 0.000199755, No Iterations 2
time step continuity errors : sum local = 4.25689e-07, global = -2.24356e-09, cumulative = -2.25318e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0427183, Final residual = 0.000117493, No Iterations 2
time step continuity errors : sum local = 2.50063e-07, global = -4.87422e-09, cumulative = -2.30193e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0427806, Final residual = 0.000118101, No Iterations 2
time step continuity errors : sum local = 2.51322e-07, global = -4.84618e-09, cumulative = -2.35039e-07

(...)

PIMPLE: iteration 198
smoothSolver:  Solving for alpha.water, Initial residual = 0.00473282, Final residual = 8.68777e-10, No Iterations 3
Phase-1 volume fraction = 0.998807  Min(alpha.water) = 1.54224e-07  Max(alpha.water) = 1.0002
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.998807  Min(alpha.water) = -0.000778275  Max(alpha.water) = 1.0002
GAMG:  Solving for p_rgh, Initial residual = 0.0242179, Final residual = 1.44059e-05, No Iterations 2
time step continuity errors : sum local = 3.02184e-08, global = -2.22393e-10, cumulative = -1.3677e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0168498, Final residual = 1.09407e-05, No Iterations 2
time step continuity errors : sum local = 2.29456e-08, global = -2.17156e-10, cumulative = -1.36987e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0168702, Final residual = 1.08204e-05, No Iterations 2
time step continuity errors : sum local = 2.26937e-08, global = -2.14002e-10, cumulative = -1.37201e-07
PIMPLE: iteration 199
smoothSolver:  Solving for alpha.water, Initial residual = 0.00473242, Final residual = 1.01542e-09, No Iterations 3
Phase-1 volume fraction = 0.998807  Min(alpha.water) = 1.50747e-07  Max(alpha.water) = 1.0002
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.998807  Min(alpha.water) = -0.00100276  Max(alpha.water) = 1.0002
GAMG:  Solving for p_rgh, Initial residual = 0.0231552, Final residual = 1.27849e-05, No Iterations 2
time step continuity errors : sum local = 2.60952e-08, global = 7.02594e-11, cumulative = -1.37131e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0160877, Final residual = 9.21079e-06, No Iterations 2
time step continuity errors : sum local = 1.8809e-08, global = 7.69685e-11, cumulative = -1.37054e-07
GAMG:  Solving for p_rgh, Initial residual = 0.0161063, Final residual = 9.18537e-06, No Iterations 2
time step continuity errors : sum local = 1.87572e-08, global = 7.68037e-11, cumulative = -1.36977e-07
PIMPLE: iteration 200
smoothSolver:  Solving for alpha.water, Initial residual = 0.00473295, Final residual = 8.68842e-10, No Iterations 3
Phase-1 volume fraction = 0.998807  Min(alpha.water) = 1.54225e-07  Max(alpha.water) = 1.0002
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.998807  Min(alpha.water) = -0.0007774  Max(alpha.water) = 1.0002
GAMG:  Solving for p_rgh, Initial residual = 0.0384742, Final residual = 4.37332e-05, No Iterations 2
time step continuity errors : sum local = 3.03523e-07, global = 2.45412e-10, cumulative = -1.36732e-07
GAMG:  Solving for p_rgh, Initial residual = 0.00193055, Final residual = 7.99627e-06, No Iterations 2
time step continuity errors : sum local = 5.51445e-08, global = 6.02003e-10, cumulative = -1.3613e-07
GAMG:  Solving for p_rgh, Initial residual = 0.000402165, Final residual = 6.91129e-09, No Iterations 19
time step continuity errors : sum local = 4.76783e-11, global = 1.11926e-11, cumulative = -1.36119e-07
PIMPLE: not converged within 200 iterations

Notable is that (1) alpha is out of its allowed bounds [0,1], which did not happen without PIMPLE,
and (2) the solution does not (or barely) improve with successive outer iterations. Then, remarkably, in the very last iteration convergence is great! This is presumably because under-relaxation is disabled. This is something I have seen several times before, but I have never understood.

The relevant files are listed below.
I should add that changing relTol for p_rgh from 1e-02 to 1e-03 yields a slightly better result, which still is not close to converging after 200 outer iterations.
Reducing the under-relaxation factors from 0.3 to 0.1 is yet much better, allowing PIMPLE to converge in 130 iterations, which ... is still an awful lot of iterations. Furthermore, over time it increased to 180 iterations.
Further reducing the under-relaxation factor from 0.1 to 0.01 yields a better residual (~5e-3) after 200 iterations, but PIMPLE still did not converge.

Question:
But then, how can PIMPLE ever yield a better(=faster) result than having only one outer iteration...?
In other words, what am I doing incorrectly?

Relevant Files:
controlDict:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     interFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.070;

deltaT          1e-6;

writeControl    adjustableRunTime;

writeInterval   5e-4;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

adjustTimeStep  yes;

maxCo           2;
maxAlphaCo      2;
maxDeltaT       1;

libs
(
        "libsimpleFunctionObjects.so"

        "libsimpleSwakFunctionObjects.so"
        "libswakLagrangianParser.so"
        "libswakFunctionObjects.so"
        "libsimpleLagrangianFunctionObjects.so"
);

functions
{
    #include "bubPositionDict" //Monitors the position of the bubble over time
}

// ************************************************************************* //

fvSolution:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      2;
        nAlphaSubCycles 1;
        cAlpha          1;

        MULESCorr       yes;
        nLimiterIter    2;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-07;
        relTol          0;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          1e-02;
        smoother        DIC;
        //nPreSweeps      0;
        //nPostSweeps     2;
        //nFinestSweeps   2;
        cacheAgglomeration no;
        //nCellsInCoarsestLevel 10;
        //agglomerator    faceAreaPair;
        //mergeLevels     1;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    "pcorr.*"
    {
        $p_rghFinal;
        tolerance       1e-04;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-07; // Doesn't do anything with momentumPredictor=no
        relTol          1e-03; // Doesn't do anything with momentumPredictor=no
        nSweeps         1;
    }

    UFinal
    {
        $U;
        relTol          0;
    }

    "(k|B|nuTilda)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-08;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor no;
    nOuterCorrectors 200;
    nCorrectors     3;
    nNonOrthogonalCorrectors 0;

    pRefPoint      (0 0 0);
    pRefValue      0;

    residualControl
    {
        U
        {
                tolerance  1e-06;
                relTol     0;
        }
        p_rgh
        {
                tolerance  1e-07;
                relTol     0;
        }
     }
}

relaxationFactors {
   fields
   {
       p_rgh        0.3;
       p_rghFinal   1;
   }
   equations
   {
       "alpha.*|U|k|epsilon"            0.3;
       "(alpha.*|U|k|epsilon)Final"     1;
   }
}


// ************************************************************************* //

fvSchemes:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    div(rhoPhi,U)   Gauss upwind;
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss linear;
//    div((muEff*dev(T(grad(U))))) Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

// ************************************************************************* //

-------

Thanks in advance.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
PIMPLE – the value of the final under-relaxation factor	Zbynek	OpenFOAM	9	December 22, 2023 05:26
pimpleDyMFoam computation randomly stops	babapeti	OpenFOAM Running, Solving & CFD	5	January 24, 2018 05:28
error while running modified pimple solver	R_21	OpenFOAM Programming & Development	0	May 28, 2015 06:59
A question on the PIMPLE algorithm	GerhardHolzinger	OpenFOAM Running, Solving & CFD	4	February 13, 2015 06:49
Help for the small implementation in turbulence model	shipman	OpenFOAM Programming & Development	25	March 19, 2014 10:08