[FloB]

Hi everyone,

I want to simulate the propagation of a Stokes V wave (T=7s, H=6m, depth=29m) thanks to olaIsoFlow solver. My mesh is very simple : it is a 2D rectangle 500 m * 60 m, I'm quite sure it is correct.

When I use PISO algorithm (nOuterCorrector=1), results are quite close to theory : sea elevation and velocity profiles are consistent.
I decided to try PIMPLE to see if results could be better, I read many threads on the forum plus this article. I set nOuterCorrectors to 50.

Unfortunately, my results are very inconsistent : I tried several relaxation factors for pressure and velocity (0.05, 0.1, 0.2, 0.4, 0.6, 0.8 and 1 for both) and despite an improvement for small relaxation factors, it is not sufficient as you can on the image attached. Residual tolerances (that is to say tolerances of residualControl in the fvSolution file) are set to 5e-3 for both variables.


I realized that the smaller was the relaxation coefficient, the shorter the calculation lasted. I don't understand it because I read that under relaxation led to slow convergence == longer calculation ?


I also tried to change residual tolerance (I put always the same tolerance for pressure and velocity) : 1e-4 and 5e-4 were apparently too ambitious (PIMPLE did not converge after 50 iterations (=nOuterCorrector) ).
However, convergence could be reach if I decrease maxCo number, but as it is already small (0.35) I would prefer to find another solution, especially that it still provides inconsistent results.

For tolerance=1e-3, it converged after 2 or 3 iterations for relax factors < 0.4, but for relax factors > 0.4, it does not converge for any time step, leading to a long calculation time.
For tolerance=5e-3 it always converged in 2 iterations, exept for few cases with relax factor = 1.


So it seems that most of the time, two extrem cases happen :
1) At the first iteration Ux, Uz and p_rgh have their initial residuals smaller than residual tolerances, so there will be an iteration 2 which will be the final one.
2) At the first iteration, initial residuals are too big, and along PIMPLE iterations, initial residuals do not change sufficiently to be inferior to residual tolerances : PIMPLE does not converge.


As we saw, for small relaxation factors I am in case 1) : but why results are so different from theory when there are 2 iterations, whereas when there is only 1 iteration (nOuterCorrector=1) results are ok ?

I also tried, to fix small value of nOuterCorrector (2, 3, 5...) without any residual tolerance : it is always inconsistent. I changed nCorrectors, it did not change any of my conclusions.


Below, are some files corresponding to relaxation factor = 0.2 and residual tolerance = 5e-3

fvSolution

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

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      1;
        nAlphaSubCycles 3;
        cAlpha          1;
        alphaOuterCorrectors  yes; 



        isoFaceTol      1e-10;
        surfCellTol     1e-6;
        nAlphaBounds    3;
        snapTol         1e-12;
        clip            true;
        writeSurfCells  false;
        writeBoundedCells false;
        writeIsoFaces   false;



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

    pcorr
    {
        solver          GAMG;   
        preconditioner  DIC;
        tolerance       1e-6;
        relTol          0.01;
	smoother	DIC;
    }

    pcorrFinal
    {
        relTol  0;
        preconditioner  DIC;
        smoother  DIC;
        tolerance  1e-6;
        solver  GAMG;     
    }

    p_rgh
    {
        solver          GAMG;   
        tolerance       1e-06;
        relTol          0.01;
	smoother	DIC;
        cacheAgglomeration no;        
    }

    p_rghFinal
    {
        $p_rgh;
        relTol  0;
        tolerance  1e-6;
    }


    U
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
    }

    UFinal
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0;
    }





}

PIMPLE
{
    momentumPredictor   on;
    nOuterCorrectors    50;  
    nCorrectors         1;
    nNonOrthogonalCorrectors 1;   
    pRefCell  0;
    pRefValue  0;

    residualControl
    {
	U
	{
		tolerance	5e-3;
		relTol		0;
	}
	p_rgh
	{
		tolerance	5e-3;
		relTol		0;
	}
    }

}

relaxationFactors
{

    fields
    {
	p_rgh	0.2;
	p_rghFinal	1;
    }
    equations
    {
        U	0.2;
	UFinal	1;
    }
}



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

fvSchemes

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  plus                                  |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     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 limitedLinearV 1;
    div(U)  Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;

}

laplacianSchemes
{
    default         Gauss linear corrected; 
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;  
}


fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha.water;
}

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

ControlDict

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

application     olaIsoFlow;

startFrom       latestTime;

startTime       0.0;

stopAt          endTime;

endTime         70;

deltaT          0.0001;  

writeControl    adjustableRunTime;

writeInterval   0.2;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;  

writeCompression off;

timeFormat      general;

timePrecision   6;   

runTimeModifiable yes;

adjustTimeStep  yes;

maxCo           0.35;
maxAlphaCo      0.35;

maxDeltaT       0.1;

functions
{
//    #include "probe"

    line
    {
        type            sets;
    	libs ("libsampling.so");
        enabled         true;
      	writeControl    writeTime;
	writeInterval   1;
       	fixedLocations false;
	interpolationScheme cellPoint;
        setFormat       raw;
        sets
       	(
	probe0
	{
		type uniform;
		axis distance;
		start   ( 0.000 0.500 0.000 );
		end     ( 0.000 0.500 40 );
		nPoints 201;
	}
	probe10
	{
		type uniform;
		axis distance;
		start   ( 10.000 0.500 0.000 );
		end     ( 10.000 0.500 40 );
		nPoints 201;
	}

	probe150
	{
		type uniform;
		axis distance;
		start   ( 150.000 0.500 0.000 );
		end     ( 150.000 0.500 40 );
		nPoints 201;
	}

	probe250
	{
		type uniform;
		axis distance;
		start   ( 250 0.500 0.000 );
		end     ( 250.000 0.500 40 );
		nPoints 201;
        }
	probe290
	{
		type uniform;
		axis distance;
		start   ( 290 0.500 0.000 );
		end     ( 290.000 0.500 40 );
		nPoints 201;
        }
	probe340
	{
		type uniform;
		axis distance;
		start   ( 340 0.500 0.000 );
		end     ( 340 0.500 40 );
		nPoints 201;
        }
	probe400
	{
		type uniform;
		axis distance;
		start   ( 400 0.500 0.000 );
		end     ( 400.000 0.500 40 );
		nPoints 201;
        }
	probe490
	{
		type uniform;
		axis distance;
		start   ( 490 0.500 0.000 );
		end     ( 490.000 0.500 40 );
		nPoints 201;
        }
        );
        fields
        (
            alpha.water
	    U
        );
    }

};


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

Can someone help me to understand why PIMPLE does not work for me ?

Thanks in advance,

Florian

[FloB]

I tried to use olaFlow (which uses MULES for VOF resolution) instead of olaIsoFlow (which uses isoAdvector) and results are much better.

It is consistent because this time, running time increases when I decrease relaxation factor.

In addition, I can decrease residual tolerance to smaller values like 1e-5 and keep PIMPLE convergence (under 50 outer iterations).


Despite the fact that now I can reach really small initial residuals, results are again less closed to theory than when I choose nOuterCorrectors=1 (with isoAdvector or MULES).
Indeed, there is phase shift, which increases in time, between theory and my simulation, plus an important wave height dissipation. These problems are not created at the inlet, they appear along the propagation.


So my questions are :

Does anyone know why using PIMPLE with isoAdvector does not work for me whereas it works (even if it is not optimal) with MULES ?

Why using nOuterCorrector>1 create a phase shift (I tried several couples of relaxation factors) with theory whereas it is not the case for nOuterCorrectors=1 ?



Thanks in advance,

Florian

[pferro]

Dear Florian.


Complex topic because many settings play a key role : mesh density, schemes, co number, time step..


This is very well detailed in this reference paper : https://arxiv.org/pdf/1804.01158.pdf.
Note that they also obtain bad result with nOuterCorrector > 1 for some settings.



I have done the same kind of tests on a periodic domain and I obtained very similar results. I was able to reproduce the bad behavior when running in PIMPLE mode.


Maybe some ideas :

- try to use a 2nd order schemes for ddt
- try to switch to linearUpwindV grad(U) for div(rhoPhi,U) and add a limiter on the velocity gradient : grad(U) cellLimited Gauss linear 1. Try also without limiter.

-you are using relaxation factor for p_rgh. Well I am rather septic for that. From experience I never relax p_rgh when using interFoam-like solvers.
-your tolerance on p_rgh solver is too high. use 1e-8.
-use a very low Co number : in the paper they use (if I remember) 0.05...
-don't be to focus on residuals.


Good luck


Paulin

[HPE]

I am not sure PIMPLE could provide better level of accuracy or stability for a given problem which seems to be a difficult numerical one in comparison to, say, single-phase incompressible flows.

PISO is mathematically consistent, published in a peer-reviewed journal. PIMPLE is a heuristic method to provide a faster yet stable computations for industrial customers. No mathematical proof no peer-reviewed journal paper. Yet very useful and powerful at times.

If you are coming from an academic habitat, I strongly recommend you to stick with PISO.

These are just my speculative opinions, however. Don't have any evidence. Just experience and gut feeling (and strong support for peer reviews).

Hope these help a bit.

[FloB]

Dear Paulin,

Thank you for the paper but I have already read it : it is really interesting indeed. I made the same study as the one mentioned in the paper but only with PISO (nOuterCorrectors=1).

Therefore, maybe I could improve my results if I try to do the same with PIMPLE but I doubt.

I'll try your suggestions, but concerning maxCo, it would be disappointing if I had to work below 0.1 : if I understood well, normally PIMPLE is made to earn time by working with higher maxCo (even >1). Moreover, my PIMPLE simulations can last three or four times more than my PISO simulations with maxCo=0.35 : imagine with maxCo=0.05, we completly loose the interest of PIMPLE.
__________________________________________________ ___________

Dear HPE,

Thank you for your insight, I agree with you, it seems that working with PISO is the best solution for me.

[pferro]

Dear FloB,

Indeed, if you have good result with PISO mode you can work with it.

It's true that PIMPLE algorithm allows you to have larger time step.
But considering wave propagation simulations maybe you can think more in terms of time step per wave period. Maybe you could work with fixed time step and perform a sensitivity analysis on it. Let's say something like 400, 200 and 100 time steps per wave period.

As you can see there are plenty of ways for working on this topic.

Paulin