Difference between revisions of "OpenFOAM guide/The PIMPLE algorithm in OpenFOAM"

From OpenFOAMWiki
(Added new reference)
(Changed the description - more information - better english and more details)
Line 1: Line 1:
 
== New reference ==  
 
== New reference ==  
If you want to check out the latest work from Tobias Holzmann about the usage of the PIMPLE algorithm, you can checkout his book »Mathematics, Numerics, Derivations and OpenFOAM(R)«. Here, you will find an own chapter about the whole stuff. In addition reference [21] on www.holzmann-cfd.de explain the pressure-velocity coupling very well.
+
If you want to check out the latest work from Tobias Holzmann about the usage of the PIMPLE algorithm, you can checkout his book »Mathematics, Numerics, Derivations and OpenFOAM(R)«. In the book you will find an own chapter which explains the algorithm in detail. In addition reference [21] on www.holzmann-cfd.de explains the pressure-velocity coupling very well. The following explanation is an old blog entry from Tobias Holzmann which is not anymore kept up to date.
  
== Storyboard ==
+
== The test case ==
First I wanted to write a nice story like »Ignaz« but that would only copy the good attention of the author and of course I am not a good story-teller. Therefore I will go directly to the topic and will tell you how the PIMPLE algorithm is working in OpenFOAM, how to set up the keywords and which abilities you have to control the algorithm. The test case it the pitzDaily case in the pimpleFoam folder of incompressible fluids. Only the Courant number was decreased to 3.
+
The test case it the pitzDaily tutorial that comes with the standard OpenFOAM installation/compilation. The case can be found in the pimpleFoam folder of incompressible fluids. The case was modified in order to allow a maximum Courant number of 3. The whole procedure was done with OpenFOAM-2.3 and is not anymore up to date.
  
  
 
== Case A ==
 
== Case A ==
In the first case we are using the tutorial case and change the fvSolution file. In the pimple control directory we will remove all entries to get only a dummy body.
+
At the beginning, we are using the tutorial case and change the fvSolution file. Within the pimple control directory we remove all entries. The outcome is a so called dummy:
 
<bash>
 
<bash>
 
PIMPLE
 
PIMPLE
Line 13: Line 13:
 
}
 
}
 
</bash>
 
</bash>
After starting the simulation you can see that the solver is still running but after a while it will crash. The following questions occur at that case:
+
After starting the simulation we realize that the solver does not crash and runs. However, based on the change of  Co, the solver will crash after a while. The following questions arise now:
 
#why is the solver working without complaining that keywords are missing
 
#why is the solver working without complaining that keywords are missing
 
#in which mode is the solver running when there are no keywords
 
#in which mode is the solver running when there are no keywords
Line 20: Line 20:
 
:'''Question 1''':
 
:'''Question 1''':
  
:If you don't insert any keywords into the PIMPLE control dictionary, the solver is still running using the default settings. After calling the pimpleControl the solver initialize the necessary keywords with default values. These can be checked out in the ''pimpleControl.C'' file:
+
:If we don't insert any keyword to the PIMPLE control dictionary, the solver is using the default value settings. To understand this, we have to go to the source code and being able to analyze the C++ code. While starting the application (pimpleFoam) we generate a new object of the pimpleControl class. Analyzing this peace of source code, we can see that there are pre-defined (so called lookupOrDefault) values.  To demonstrate this, a snippet of the ''pimpleControl.C'' file is given below:
 
  <cpp>
 
  <cpp>
 
   121 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //  
 
   121 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //  
Line 35: Line 35:
 
   132    read();
 
   132    read();
 
</cpp>
 
</cpp>
:As you can see there are several keywords that can be used if you are using the PIMPLE algorithm. These keywords are set »by default« to the values you can see at the code above. After initialize the values, the constructor will call the »read()« function which will check out the user defined values and overwrite the above values with user values or other values which are necessary that the solver is working. If there are no values in the PIMPLE control available the solver is using the default one. The functions is given below:
+
:As we can see, there are several keywords that are initialized with default values. After the initializing phase, the constructor is calling the »read()« function which checks if the user defined the necessary keywords to overwrite the default values. If no keyword is set, the default values are set (see listing below):
 
<cpp>
 
<cpp>
 
     39 void Foam::pimpleControl::read()     
 
     39 void Foam::pimpleControl::read()     
Line 49: Line 49:
 
     49 }
 
     49 }
 
</cpp>
 
</cpp>
:As you can see here, this function is looking in the '''fvSolutions''' file for some entries in the PIMPLE control structure:
+
:As we can see, the read function is looking in the '''fvSolutions''' file for the corresponding entries in the PIMPLE control dictionary:
 
<pre>
 
<pre>
 
nOuterCorrctors (nCorrPIMPLE) is set by default to 1
 
nOuterCorrctors (nCorrPIMPLE) is set by default to 1
Line 61: Line 61:
 
:'''Question 2''':
 
:'''Question 2''':
  
:If you are running the PIMPLE solver without any keywords it is operating in PISO mode and therefore you need very small time steps. Thats the reason why the simulation will blow up with these settings (Co = 3). That the solver is operating in PISO mode can be shown in the source code too, and will report it on the screen:
+
:If we are running the PIMPLE solver without any keywords, the solver operatings in the so called PISO mode. Therefore, we need to fulfill the stability criterion Co < 1. Based on the fact that we set CoMax to 3, this is  the reason why the simulation will blow up. OpenFOAM outputs the operating mode of the solver after starting the solver:
 
<cpp>
 
<cpp>
 
   134    if (nCorrPIMPLE_ > 1)   
 
   134    if (nCorrPIMPLE_ > 1)   
Line 77: Line 77:
  
  
:'''Residuals''' if you are using a empty PIMPLE directory
+
:'''Residuals:''' If we use an empty PIMPLE directory the residual plot looks as follow:
  
 
:[[Image:Case 1 Residual.png]]
 
:[[Image:Case 1 Residual.png]]
:As you can see the residuals are not converging and the solution blow up. The reason is due to the high Co number (Co = 3) and the fact that without keywords the PIMPLE algorithm is working in PISO mode. To solve the problem you have to reduce the Co number but that is not is not the aim to introduce merged PISO - SIMPLE algorithm.
+
:As we see, the residuals do not converge and the solution blow up after round about 27 time steps (each iteration = one time step). Again, the blow up is based on the limitation of the PISO algorithm which forces us to keep Co less than 1.  There are two options to get the simulation stable. The first one is to limit the Courant number, which is in fact not the topic we are talking about here. The second one is to use the PIMPLE algorithm in the correct way. The PIMPLE algorithm combines the PISO with the SIMPLE. That means, the problem is transient but we use the SIMPLE (steady-state) treatment to find the steady-state solution for each time step. This allows us to increase the Courant number larger than one.
  
 
:'''Legend''':
 
:'''Legend''':
  
::* Iteration at the bottom means one PISO iteration; that means in each iteration you go on with the time. In other words every iteration solves a new time step.
+
::* Iteration at the bottom means one PISO iteration; that means in each iteration we move on by delta t. In other words every iteration solves a new time step.
  
  
 
== Case B ==
 
== Case B ==
In the second case we are using the standard entries in the '''fvSolutions''' file used in the tutorial.
+
The second case sets the first keyword in order to change the operating mode of the solver. We change the dummy directory as follow:
 
<cpp>
 
<cpp>
 
PIMPLE
 
PIMPLE
Line 95: Line 95:
 
}
 
}
 
</cpp>
 
</cpp>
Now we correct the pressure field twice and that means that we calculate the following equation twice:
+
Now we correct the pressure field twice which means that the following equation is calculated twice:
  
 
<math>
 
<math>
Line 101: Line 101:
 
</math>
 
</math>
  
In that case the pressure - momentum coupling is stable enough to work and the solver will not blow up.
+
In our case B,  the pressure - momentum coupling is now stable enough to produce a non-diverging solution. Thus, the solver is not blowing up.
  
 
:'''Advantage'''
 
:'''Advantage'''
::* While you are using nCorrectors with a number higher than 1 the PIMPLE algorithm is still working in PISO mode with the advantage that you calculate your pressure twice, this is also available in the normal PISO solver but is still not the aim of the PIMPLE solver.
+
::* While using nCorrectors with a number higher than one, the application (pimpleFoam) is still working in PISO mode with the advantage that we re-calculate the pressure twice. This PISO loop option is always available in all transient FOAM solvers, if the application uses the PISO or PIMPLE control class.
  
  
:'''Residuals''' if you are using the PIMPLE algorithm with {{tt|nCorrectors 2}}
+
:'''Residuals:''' If we have a look at the new residual plots in which we have {{tt|nCorrectors 2}}, we get the following:
 
:[[Image:Case 2 Residual.png]]
 
:[[Image:Case 2 Residual.png]]
  
:The residual of the pressure is much better and the solution is stable in that moment - not tested what happen if you solve more seconds.
+
:The residuals of the pressure are much better than in case A and the solution is stable for the shown iterations - not tested if the solver will actually blow up, if we run the case for longer time.
  
 
:'''Legend:'''
 
:'''Legend:'''
::* For the momentum one iteration means »one time step«, in other words that means in that case we were solving round about 65 time steps
+
::* For the momentum lines, each iteration represents »one time step«. Thus, we are solving around 65 time steps
::* For the pressure we solve the equation given above two times and therefore it is not anymore possible to say that one iteration is one time step. In that case two iterations are necessary to go on in time.
+
::* For the pressure we solve the equation given above twice and therefore, it is not anymore possible to say that one iteration is one time step. Of course, it is obvious that two iterations represents one time step but it is not always the case.
  
  
 
==  Case C ==
 
==  Case C ==
 +
Now we are going to introduce the real merged PISO - SIMPLE algorithm and therefore we are using the advantage of the SIMPLE mode for each time step (relaxation). What does that mean? If we are using the PISO mode, we are limited to the time step which  fulfills the condition that Co < 1. Based on the definition of the Courant number one knows that this can lead to very expensive computational resources while solving a real time problem. To speedup the simulation we have to overcome the limit of Co < 1 which can be achieved by using the SIMPLE algorithm. The reason for that is based on the fact that the SIMPLE algorithm is developed for the final, steady-state, condition. To reach that flow behavior rapidly, we are actually not interested in the transient behavior and development of the flow (ddt terms = zero) and thus, dt should be as large as possible. Of course there is much more behind which is not mentioned right now. More details can be found in the book mentioned at the beginning or in other literature.
  
Now we will use the MERGED PISO - SIMPLE algorithm and therefore we are using the advantage of the SIMPLE mode idea (relaxation). What does that mean? In fact it is really easy. If you are using the PISO mode you are strict limited to the time step (Co < 1). That mean that it is very expensive to solve a real time problem in transient manner for long time especially in very complex geometries. SIMPLE is developed for steady-state condition to reach that state very fast but (normally) does not contain time derivation and therefore no time information. To provide a stable simulation you have to introduce relaxation factors that limit the values of the new values.  
+
However, a very large time step or actually neglecting the temporal derivative changes the characteristic of the equations. For such cases the SIMPLE algorithm families (SIMPLE, SIMPLEC, SIMPLER...) were developed. Based on the fact that in the normal SIMPLE algorithm, a pressure term is neglected, we need to have under-relaxation (even with SIMPLC it is needed in most cases but the relaxation factors can be chosen larger). In addition it should be mentioned that there are two different relaxation methods available in OpenFOAM: Field and matrix relaxation. More about that is given in the book mentioned at the beginning.
 
+
 
+
<math>
+
\phi^n = \phi^{n-1} + \alpha_\phi ( \phi^{\mathrm{new}} - \phi^{n-1})
+
</math>
+
 
+
 
+
In other words, the new value <math>\phi^n</math> is a combination of the old value <math>\phi^{n-1}</math> and the new one <math>\phi^{\mathrm{new}}</math> which is calculated. This make the simulation stable but not time conservative. It allows you to reach the steady-state condition very fast; time information is not needed.
+
 
+
  
 
:'''Miss understanding'''
 
:'''Miss understanding'''
  
:A lot of people miss understand the function of the PIMPLE algorithm. The problem in the case we are using here is, that this case is very simple. Imagine a very complex simulation in which it is very straight forward to simulate small time steps. If you want to have bigger time steps you have to increase your Co number which leads you to unstable convergence. At that moment you are in a discrepancy because you want go on faster in time step to reduce computational effort but it is not possible because you are limited in time step due to the Co number.
+
:A lot of people miss understand the capability and usage of the PIMPLE algorithm. The main problem in the case we are analyzing right now is, that this is very simple and does not need special numerical treatments. However, in complex geometries and with different models which causes stiff problems, the PIMPLE algorithm is the method of choice. Imagine a very complex CFD simulation with small cell sizes and high velocities, it is straight forward to use small time steps in the PISO mode. However, If one wants to have bigger time steps one has to increase the Co number which break the criterion of Co < 1. At that moment we get a discrepancy because one one hand we want go on faster in time in order to reduce the computational effort but it is not possible because we are limited based on C o < 1. Thats why we use the PIMPLE mode then. Keep in mind that each time step can be much more expensive than a time step in the PISO calculation but the time step itself can be much larger. In addition it can occur that we loose significant information of the flow which changes the flow pattern.
  
  
 
=== Case C.1 ===
 
=== Case C.1 ===
Now we will use the following '''fvSolutions''' file:
+
For the C1 case we change the '''fvSolutions''' file as follow:
 
<cpp>
 
<cpp>
 
PIMPLE
 
PIMPLE
Line 144: Line 136:
 
}
 
}
 
</cpp>
 
</cpp>
Compared to Case B there is not a big difference. The only difference is that we calculate the pressure - momentum coupling in one time step twice:
+
Compared to case B, we added the nOuterCorrectors keyword and set to two while changing the pressure correction back to one again. The new keyword and the value > 1 activates the PIMPLE loop and lead to a re-calculation of the pressure - momentum coupling (here two times) in one time step.
:# Calculate momentum
+
:# Construct momentum matrix
 +
:# Construct the pressure matrix using the momentum matrix
 
:# Calculate pressure
 
:# Calculate pressure
:# Recalculate momentum from new pressure
+
:# Correct the velocities with the new pressure field
:# Recalculate pressure from new momentum
+
:# Reconstruct the momentum matrix with the new velocities
 +
:# Construct the pressure matrix using the momentum matrix with the updated velocities
 +
:# Calculate the pressure
 +
:# Correct the velocities with the new pressure field
  
 
::1.    <math>\frac{U}{\partial t} + \nabla \cdot (UU) + \nabla \cdot R = -\nabla p</math>
 
::1.    <math>\frac{U}{\partial t} + \nabla \cdot (UU) + \nabla \cdot R = -\nabla p</math>
Line 158: Line 154:
 
::<math>\rightarrow</math> correct U with <math>p_{\mathrm{new}} \rightarrow U_\mathrm{corrected}</math>
 
::<math>\rightarrow</math> correct U with <math>p_{\mathrm{new}} \rightarrow U_\mathrm{corrected}</math>
  
::<math>\rightarrow</math> now you use the <math>U_\mathrm{corrected}</math> and <math>p_{\mathrm{new}}</math> for the second loop
+
::<math>\rightarrow</math> now we have the <math>U_\mathrm{corrected}</math> and <math>p_{\mathrm{new}}</math> for the second loop
  
 
::3.    <math>\frac{U}{\partial t} + \nabla \cdot (UU) + \nabla \cdot R = -\nabla p</math>
 
::3.    <math>\frac{U}{\partial t} + \nabla \cdot (UU) + \nabla \cdot R = -\nabla p</math>
Line 169: Line 165:
  
  
During this calculation it has not a big advantage and is like Case A calculated twice and then move on in time step (with the only difference, that the turbulence is calculated at the end. That means in PIMPLE Iteration 2 due to the fact that ''turbOnFinalIterOnly = true''). You can assume that the calculation will blow up which is shown in the residual plot.
+
Case C1 does not has big advantages and is mostly like case A with the difference that we accelerate error accumulation and calculate the turbulence equations only at the end of the PIMPLE loop (here after two pimple iterations based on the fact that ''turbOnFinalIterOnly" is set to true as default.
  
  
:'''Residuals''' if you are using the PIMPLE algorithm with the settings above
+
:'''Residuals:''' Checking the residual plots give us:
 
:[[Image:Case 3a Residual.png]]
 
:[[Image:Case 3a Residual.png]]
:As you can see it is much worse compared to Case 1. At that moment there is none advantage of the pimple mode that we are using but we are using it in a wrong way.
+
:As already mentioned, it is amplifying the errors and the divergence occur earlier compared to case A. Keep in mind that the legend here represents the outer loops. Thus, the simulation diverges already after half of the time compared to case A. The reason is that we use the PIMPLE algorithm in a wrong way.
  
 
=== Case C.2 ===
 
=== Case C.2 ===
  
Now we want to generate a complete PIMPLE case. The '''fvSolutions''' file looks like that:
+
Case C.2 will set-up the complete PIMPLE mode. Therefore, we have to add the relaxationFators dictionary for the fields and equations. Thus, the '''fvSolutions''' file looks as follow:
 
<cpp>
 
<cpp>
 
PIMPLE
 
PIMPLE
Line 184: Line 180:
 
     nNonOrthogonalCorrectors 0;
 
     nNonOrthogonalCorrectors 0;
 
     nCorrectors          1;
 
     nCorrectors          1;
     nOuterCorrectors    2;  // Set to 50  
+
     nOuterCorrectors    50;
 
+
}
  
 
relaxationFactors
 
relaxationFactors
Line 202: Line 198:
 
</cpp>
 
</cpp>
  
We introduced the relaxation factors to obtain a smooth convergence. The problem now is that we use the relaxation factors in the first pimple loop but after we reach the last one (here iteration 2 because ''nOuterCorrectors = 2'') the relaxation factor will be 1 again to ensure time consistence but this makes no sense in the following case. Therefore we will change the nOuterCorrectors to 50. Now we calculate 50 times the pressure - momentum coupling. Therefore 49 times with relaxation (like SIMPLE) and the last time without relaxation. For cases that are very unstable this is perfect to get a smooth convergence. In other words we are using the idea of SIMPLE mode (under relaxation) between the time steps to converge the solution at these time step to a certain limit. The limit is set in a way to be sure the solution is accurate in time and space.
+
We introduced the relaxation factors for the SIMPLE calculation in between the single time steps to obtain a smooth convergence rate. The application relaxes the fields and equations for each pimple loop as given above. After we reach the 50 outer corrections, the final flag is set and the relaxation factors are set to one again. A discussion about that is given in [1].  However, in case C.we calculate the pressure - momentum coupling now for 50 times (49 times with relaxation and the last time without relaxation - if  the .*Final flag is set to 1 or neglected). For stiff problems or/and weak described systems (boundary conditions etc.) this method allows us to get a smooth convergence rate.
  
:'''Residuals''' if you are using the right PIMPLE algorithm with relaxation
+
:'''Residuals:''' Analyzing the residual plots again, we get the following:
 
:[[Image:Case 3b Residual.png]]
 
:[[Image:Case 3b Residual.png]]
:You can see that the solution is stable and the convergence is very good. Further more you can see the pimple loops; here from zero to fifty iteration, 100 to 150 and so on. After you finished the pimple loop we go on in time and do the same thing again. If you introduce the ''nCorrectors > 1'' you will calculate the pressure field n times (more hints to that later).
+
:We see that the solution is stable and the convergence rate is very good. Furthermore, we see the steady-state behavior between the single time steps; here from zero to fifty iteration, 100 to 150 and so on. After finishing the pimple loop, we move on in time and do the same calculation procedure again. In addition we could set the ''nCorrectors > 1'' which is briefly discussed/mentioned later on.
  
  
:'''Do not forget'''
+
:'''Never forget'''
  
:At the moment we have 50 Iterations per time step and therefore the computational effort could be bigger as before. To avoid that you should always add residual controls to stop the pimple loop after reaching some special criterion » Case C.3
+
:Right now we solve 50 corrections (momentum - pressure coupling) for each time step and therefore the computational effort could be much larger than before. To avoid this, we always should add the residual control in order to leave the pimple loop after we reach the residual control of each quantity we specify » case C.3
  
 
=== Case C.3 ===
 
=== Case C.3 ===
At last, we introduce the residual control to exit the pimple loop after we reached a good solution during one time step. The '''fvSolutions''' file looks like that:
+
At the end we add the residual control section to the fvSolutions file in order to leave the momentum-pressure coupling loops after everything is converged to speedup the calculation and use all advantages from the PIMPLE algorithm:
 
<cpp>
 
<cpp>
 
PIMPLE
 
PIMPLE
Line 252: Line 248:
 
</cpp>
 
</cpp>
  
:'''Residuals''' if you are using the right PIMPLE algorithm with relaxation and residualControl
+
:'''Residuals:''' A last time we are analyzing the residual plots of the calculation with relaxation and residual control:
 
:[[Image:Case 3c Residual.png]]
 
:[[Image:Case 3c Residual.png]]
:As you can see now, the time iteration is faster because after the pimple loop will reach the residual control it will go to the next time step.
+
:We see that after the first time step - which requires almost 50 outer corrections - the pimple loop is stopped earlier and therefore the simulation is faster (each time step needs less outer corrections). If the PIMPLE algorithm will converge with the settings one is using can be checked out after a new time step is started:
 
<pre>
 
<pre>
smoothSolver:  Solving for Ux, Initial residual = 2.07804e-06, Final residual = 2.07804e-06, No Iterations 0
+
.
smoothSolver:  Solving for Uy, Initial residual = 3.42321e-06, Final residual = 3.42321e-06, No Iterations 0
+
.
GAMG:  Solving for p, Initial residual = 0.000407312, Final residual = 3.14679e-08, No Iterations 8
+
.
time step continuity errors : sum local = 5.66964e-12, global = 3.09312e-13, cumulative = -8.61807e-06
+
smoothSolver:  Solving for epsilon, Initial residual = 0.00512886, Final residual = 2.08382e-06, No Iterations 2
+
 
smoothSolver:  Solving for k, Initial residual = 0.0660137, Final residual = 1.08191e-06, No Iterations 3
 
smoothSolver:  Solving for k, Initial residual = 0.0660137, Final residual = 1.08191e-06, No Iterations 3
 
'''PIMPLE: converged in 28 iterations'''
 
'''PIMPLE: converged in 28 iterations'''
 
ExecutionTime = 13.84 s  ClockTime = 14 s
 
ExecutionTime = 13.84 s  ClockTime = 14 s
 
</pre>
 
</pre>
:Although we set 50 iterations for the pimple loop it will exit the loop after the residual control settings are true. That saves computational effort.
+
:Although we set the total amount of 50 outer corrections, the calculation will exit the loop after the residual control settings are fulfilled. In the particular case it stops after 28.
  
  
 
== Some hints ==
 
== Some hints ==
  
:* the pitzDaily tutorial is to easy to show the advantage of PIMPLE in a good manner
+
:* the pitzDaily tutorial is an too easy case to show the advantage of PIMPLE in a good manner (another case is given in the book by Tobias Holzmann mentioned at the beginning)
  
:* PIMPLE is very effective in big cases to increase the Co number
+
:* The PIMPLE algorithm is very effective for stiff and large cases, for really bad initialized cases, weak boundary setup and if one wants to increase the Co number
  
:* PIMPLE speed up your simulation and could reach a stable solution due to under relaxation
+
:* The PIMPLE algorithm speeds up your simulation and allow us to reach nice convergence rates while using under relaxation in between the time steps
  
:* For high turbulent flow you should add turbOnFinalIterOnly = false to your PIMPLE control
+
:* For high turbulent flow fields we should add the keyword turbOnFinalIterOnly = false to calculate the turbulent properties each outer loop too
  
:* Set nOuterCorrectors up to a very high number (~ 50 to 200) and control your pimple loop with the residual control
+
:* Always set the nOuterCorrectors to a very high number (~ 50 to 1000) and control your pimple loop with the residual control (make sure that explizit terms converged)
  
:* Set nCorrectors to a small number (1 to 3; you have to play around)
+
:* Set nCorrectors to a small number (1 to 3)
  
:* If you set nCorrectors higher it could be possible to leave the pimple loop earlier and prevent extra calculation but could also do the opposite
+
:* If one sets nCorrectors to a high value it could be possible to leave the pimple loop earlier and prevent extra calculation but can be vice versa too
  
:* nNonOrthogonalCorr means that you correct the pressure field more often with the new calculated one:  
+
:* nNonOrthogonalCorr means that you correct the pressure field more often with the new calculated one. This is needed for high non-orthogonal meshes because we get an explizit correction term (nOuterCorrectors will do the same but with more computational effort):
  
  
Line 290: Line 284:
  
  
:* The optimum relaxation factors depend on the simulation case you are running
+
:* The optimum relaxation factors depend on the simulation case we are running and if we use the SIMPLE or SIMPLEC (for corrected) algorithm
:* The converged solution is in depended of the relaxation factors which are applied
+
:* Due to the fact that the PIMPLE algorithm allows us to use bigger time steps, it could happen that we loose important information which lead to different results
:* Due to the fact that you have bigger time steps it could happen that some physics are not resolved (extreme fast pressure drop etc.)
+
:* Right at the start of a transient simulation it is common that the application does not converge the first time steps within the given outer corrections based on the bad initial fields. This should change after a few time steps
:* At the beginning of an transient simulation it could happen that after 100 pimple loops the solution is still not converged but this will change during the simulation time
+
:* The CHT solver needs a lot of PIMPLE loops but can be reduced by tweaking the numerics (if you know how - a tutorial can be found on the website of Tobias Holzmann)
::* simple heat transfer case has at the beginning 200 pimple loops and is still not converging
+
::* after the tenth time step it converges and the pimple loops decrease significantly till it get a steady state pimple iteration of round about 3 - 5
+
  
  
'''Please feel free to change this site if I got something wrong'''
+
'''Updated 21.08.2017'''
  
 
[[Category:Numerical_methods]]
 
[[Category:Numerical_methods]]

Revision as of 23:30, 21 August 2017

1 New reference

If you want to check out the latest work from Tobias Holzmann about the usage of the PIMPLE algorithm, you can checkout his book »Mathematics, Numerics, Derivations and OpenFOAM(R)«. In the book you will find an own chapter which explains the algorithm in detail. In addition reference [21] on www.holzmann-cfd.de explains the pressure-velocity coupling very well. The following explanation is an old blog entry from Tobias Holzmann which is not anymore kept up to date.

2 The test case

The test case it the pitzDaily tutorial that comes with the standard OpenFOAM installation/compilation. The case can be found in the pimpleFoam folder of incompressible fluids. The case was modified in order to allow a maximum Courant number of 3. The whole procedure was done with OpenFOAM-2.3 and is not anymore up to date.


3 Case A

At the beginning, we are using the tutorial case and change the fvSolution file. Within the pimple control directory we remove all entries. The outcome is a so called dummy:

 
PIMPLE
{
}

After starting the simulation we realize that the solver does not crash and runs. However, based on the change of Co, the solver will crash after a while. The following questions arise now:

  1. why is the solver working without complaining that keywords are missing
  2. in which mode is the solver running when there are no keywords


Question 1:
If we don't insert any keyword to the PIMPLE control dictionary, the solver is using the default value settings. To understand this, we have to go to the source code and being able to analyze the C++ code. While starting the application (pimpleFoam) we generate a new object of the pimpleControl class. Analyzing this peace of source code, we can see that there are pre-defined (so called lookupOrDefault) values. To demonstrate this, a snippet of the pimpleControl.C file is given below:
 
   121 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * // 
   122 
   123 Foam::pimpleControl::pimpleControl(fvMesh& mesh)
   124 :
   125     solutionControl(mesh, "PIMPLE"),
   126     nCorrPIMPLE_(0),
   127     nCorrPISO_(0),
   128     corrPISO_(0),
   129     turbOnFinalIterOnly_(true),
   130     converged_(false)
   131 {   
   132     read();
As we can see, there are several keywords that are initialized with default values. After the initializing phase, the constructor is calling the »read()« function which checks if the user defined the necessary keywords to overwrite the default values. If no keyword is set, the default values are set (see listing below):
 
    39 void Foam::pimpleControl::read()    
    40 {
    41     solutionControl::read(false);
    42 
    43     // Read solution controls
    44     const dictionary& pimpleDict = dict();
    45     nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label> ("nOuterCorrectors", 1);
    46     nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
    47     turbOnFinalIterOnly_ =
    48         pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
    49 }
As we can see, the read function is looking in the fvSolutions file for the corresponding entries in the PIMPLE control dictionary:
nOuterCorrctors (nCorrPIMPLE) is set by default to 1
nCorrectors (nCorrPISO) is set by default to 1
nNonOrthogonalCorrectors (corrPISO) is set by default to 0  (not shown here)
turbOnFinalIterOnly is set by default to true
Residual control is  set to false


Question 2:
If we are running the PIMPLE solver without any keywords, the solver operatings in the so called PISO mode. Therefore, we need to fulfill the stability criterion Co < 1. Based on the fact that we set CoMax to 3, this is the reason why the simulation will blow up. OpenFOAM outputs the operating mode of the solver after starting the solver:
 
  134     if (nCorrPIMPLE_ > 1)   
  135     {
   .
   .
   .
   156     }
   157     else
   158     {
   159         Info<< nl << algorithmName_ << ": Operating solver in PISO mode" << nl
   160             << endl;
   161     }


Residuals: If we use an empty PIMPLE directory the residual plot looks as follow:
Case 1 Residual.png
As we see, the residuals do not converge and the solution blow up after round about 27 time steps (each iteration = one time step). Again, the blow up is based on the limitation of the PISO algorithm which forces us to keep Co less than 1. There are two options to get the simulation stable. The first one is to limit the Courant number, which is in fact not the topic we are talking about here. The second one is to use the PIMPLE algorithm in the correct way. The PIMPLE algorithm combines the PISO with the SIMPLE. That means, the problem is transient but we use the SIMPLE (steady-state) treatment to find the steady-state solution for each time step. This allows us to increase the Courant number larger than one.
Legend:
  • Iteration at the bottom means one PISO iteration; that means in each iteration we move on by delta t. In other words every iteration solves a new time step.


4 Case B

The second case sets the first keyword in order to change the operating mode of the solver. We change the dummy directory as follow:

 
PIMPLE
{
    nCorrectors         2;
}

Now we correct the pressure field twice which means that the following equation is calculated twice:


\nabla^2 p = f(U, \nabla p)

In our case B, the pressure - momentum coupling is now stable enough to produce a non-diverging solution. Thus, the solver is not blowing up.

Advantage
  • While using nCorrectors with a number higher than one, the application (pimpleFoam) is still working in PISO mode with the advantage that we re-calculate the pressure twice. This PISO loop option is always available in all transient FOAM solvers, if the application uses the PISO or PIMPLE control class.


Residuals: If we have a look at the new residual plots in which we have nCorrectors 2, we get the following:
Case 2 Residual.png
The residuals of the pressure are much better than in case A and the solution is stable for the shown iterations - not tested if the solver will actually blow up, if we run the case for longer time.
Legend:
  • For the momentum lines, each iteration represents »one time step«. Thus, we are solving around 65 time steps
  • For the pressure we solve the equation given above twice and therefore, it is not anymore possible to say that one iteration is one time step. Of course, it is obvious that two iterations represents one time step but it is not always the case.


5 Case C

Now we are going to introduce the real merged PISO - SIMPLE algorithm and therefore we are using the advantage of the SIMPLE mode for each time step (relaxation). What does that mean? If we are using the PISO mode, we are limited to the time step which fulfills the condition that Co < 1. Based on the definition of the Courant number one knows that this can lead to very expensive computational resources while solving a real time problem. To speedup the simulation we have to overcome the limit of Co < 1 which can be achieved by using the SIMPLE algorithm. The reason for that is based on the fact that the SIMPLE algorithm is developed for the final, steady-state, condition. To reach that flow behavior rapidly, we are actually not interested in the transient behavior and development of the flow (ddt terms = zero) and thus, dt should be as large as possible. Of course there is much more behind which is not mentioned right now. More details can be found in the book mentioned at the beginning or in other literature.

However, a very large time step or actually neglecting the temporal derivative changes the characteristic of the equations. For such cases the SIMPLE algorithm families (SIMPLE, SIMPLEC, SIMPLER...) were developed. Based on the fact that in the normal SIMPLE algorithm, a pressure term is neglected, we need to have under-relaxation (even with SIMPLC it is needed in most cases but the relaxation factors can be chosen larger). In addition it should be mentioned that there are two different relaxation methods available in OpenFOAM: Field and matrix relaxation. More about that is given in the book mentioned at the beginning.

Miss understanding
A lot of people miss understand the capability and usage of the PIMPLE algorithm. The main problem in the case we are analyzing right now is, that this is very simple and does not need special numerical treatments. However, in complex geometries and with different models which causes stiff problems, the PIMPLE algorithm is the method of choice. Imagine a very complex CFD simulation with small cell sizes and high velocities, it is straight forward to use small time steps in the PISO mode. However, If one wants to have bigger time steps one has to increase the Co number which break the criterion of Co < 1. At that moment we get a discrepancy because one one hand we want go on faster in time in order to reduce the computational effort but it is not possible because we are limited based on C o < 1. Thats why we use the PIMPLE mode then. Keep in mind that each time step can be much more expensive than a time step in the PISO calculation but the time step itself can be much larger. In addition it can occur that we loose significant information of the flow which changes the flow pattern.


5.1 Case C.1

For the C1 case we change the fvSolutions file as follow:

 
PIMPLE
{
    nCorrectors 1;
    nOuterCorrectors 2;
}

Compared to case B, we added the nOuterCorrectors keyword and set to two while changing the pressure correction back to one again. The new keyword and the value > 1 activates the PIMPLE loop and lead to a re-calculation of the pressure - momentum coupling (here two times) in one time step.

  1. Construct momentum matrix
  2. Construct the pressure matrix using the momentum matrix
  3. Calculate pressure
  4. Correct the velocities with the new pressure field
  5. Reconstruct the momentum matrix with the new velocities
  6. Construct the pressure matrix using the momentum matrix with the updated velocities
  7. Calculate the pressure
  8. Correct the velocities with the new pressure field
1. \frac{U}{\partial t} + \nabla \cdot (UU) + \nabla \cdot R = -\nabla p
2. \nabla^2 p = f(U, \nabla p)
\rightarrow p_{\mathrm{new}}
\rightarrow correct U with p_{\mathrm{new}} \rightarrow U_\mathrm{corrected}
\rightarrow now we have the U_\mathrm{corrected} and p_{\mathrm{new}} for the second loop
3. \frac{U}{\partial t} + \nabla \cdot (UU) + \nabla \cdot R = -\nabla p
4. \nabla^2 p = f(U, \nabla p)
\rightarrow p_{\mathrm{new}}
\rightarrow correct U with p_{\mathrm{new}} \rightarrow U_\mathrm{corrected}


Case C1 does not has big advantages and is mostly like case A with the difference that we accelerate error accumulation and calculate the turbulence equations only at the end of the PIMPLE loop (here after two pimple iterations based on the fact that turbOnFinalIterOnly" is set to true as default.


Residuals: Checking the residual plots give us:
Case 3a Residual.png
As already mentioned, it is amplifying the errors and the divergence occur earlier compared to case A. Keep in mind that the legend here represents the outer loops. Thus, the simulation diverges already after half of the time compared to case A. The reason is that we use the PIMPLE algorithm in a wrong way.

5.2 Case C.2

Case C.2 will set-up the complete PIMPLE mode. Therefore, we have to add the relaxationFators dictionary for the fields and equations. Thus, the fvSolutions file looks as follow:

 
PIMPLE
{
    nNonOrthogonalCorrectors 0;
    nCorrectors          1;
    nOuterCorrectors    50; 
}
 
relaxationFactors
{
    fields
    {
        p      0.3;
        pFinal   1;
    }
    equations
    {
        "U|k|epsilon"     0.3;
        "(U|k|epsilon)Final"   1;
    }
}

We introduced the relaxation factors for the SIMPLE calculation in between the single time steps to obtain a smooth convergence rate. The application relaxes the fields and equations for each pimple loop as given above. After we reach the 50 outer corrections, the final flag is set and the relaxation factors are set to one again. A discussion about that is given in [1]. However, in case C.2 we calculate the pressure - momentum coupling now for 50 times (49 times with relaxation and the last time without relaxation - if the .*Final flag is set to 1 or neglected). For stiff problems or/and weak described systems (boundary conditions etc.) this method allows us to get a smooth convergence rate.

Residuals: Analyzing the residual plots again, we get the following:
Case 3b Residual.png
We see that the solution is stable and the convergence rate is very good. Furthermore, we see the steady-state behavior between the single time steps; here from zero to fifty iteration, 100 to 150 and so on. After finishing the pimple loop, we move on in time and do the same calculation procedure again. In addition we could set the nCorrectors > 1 which is briefly discussed/mentioned later on.


Never forget
Right now we solve 50 corrections (momentum - pressure coupling) for each time step and therefore the computational effort could be much larger than before. To avoid this, we always should add the residual control in order to leave the pimple loop after we reach the residual control of each quantity we specify » case C.3

5.3 Case C.3

At the end we add the residual control section to the fvSolutions file in order to leave the momentum-pressure coupling loops after everything is converged to speedup the calculation and use all advantages from the PIMPLE algorithm:

 
PIMPLE
{
    nNonOrthogonalCorrectors 0;
    nCorrectors          1;
    nOuterCorrectors    50;
 
    residualControl
    {
        U
        {
                tolerance  1e-5;
                relTol      0;
        }
        p
        {
                tolerance  5e-4;
                relTol      0;
        }
     }
}
 
relaxationFactors
{
    fields
    {
        p      0.3;
        pFinal   1;
    }
    equations
    {
        "U|k|epsilon"     0.3;
        "(U|k|epsilon)Final"   1;
    }
}
Residuals: A last time we are analyzing the residual plots of the calculation with relaxation and residual control:
Case 3c Residual.png
We see that after the first time step - which requires almost 50 outer corrections - the pimple loop is stopped earlier and therefore the simulation is faster (each time step needs less outer corrections). If the PIMPLE algorithm will converge with the settings one is using can be checked out after a new time step is started:
.
.
.
smoothSolver:  Solving for k, Initial residual = 0.0660137, Final residual = 1.08191e-06, No Iterations 3
'''PIMPLE: converged in 28 iterations'''
ExecutionTime = 13.84 s  ClockTime = 14 s
Although we set the total amount of 50 outer corrections, the calculation will exit the loop after the residual control settings are fulfilled. In the particular case it stops after 28.


6 Some hints

  • the pitzDaily tutorial is an too easy case to show the advantage of PIMPLE in a good manner (another case is given in the book by Tobias Holzmann mentioned at the beginning)
  • The PIMPLE algorithm is very effective for stiff and large cases, for really bad initialized cases, weak boundary setup and if one wants to increase the Co number
  • The PIMPLE algorithm speeds up your simulation and allow us to reach nice convergence rates while using under relaxation in between the time steps
  • For high turbulent flow fields we should add the keyword turbOnFinalIterOnly = false to calculate the turbulent properties each outer loop too
  • Always set the nOuterCorrectors to a very high number (~ 50 to 1000) and control your pimple loop with the residual control (make sure that explizit terms converged)
  • Set nCorrectors to a small number (1 to 3)
  • If one sets nCorrectors to a high value it could be possible to leave the pimple loop earlier and prevent extra calculation but can be vice versa too
  • nNonOrthogonalCorr means that you correct the pressure field more often with the new calculated one. This is needed for high non-orthogonal meshes because we get an explizit correction term (nOuterCorrectors will do the same but with more computational effort):


\nabla^2 p = f(U, \nabla p) \rightarrow p_{\mathrm{new}} \rightarrow \nabla^2 p_{\mathrm{new}} = f(U, \nabla p_{\mathrm{new}} ) \rightarrow till nNonOrtogonalCorr is reached


  • The optimum relaxation factors depend on the simulation case we are running and if we use the SIMPLE or SIMPLEC (for corrected) algorithm
  • Due to the fact that the PIMPLE algorithm allows us to use bigger time steps, it could happen that we loose important information which lead to different results
  • Right at the start of a transient simulation it is common that the application does not converge the first time steps within the given outer corrections based on the bad initial fields. This should change after a few time steps
  • The CHT solver needs a lot of PIMPLE loops but can be reduced by tweaking the numerics (if you know how - a tutorial can be found on the website of Tobias Holzmann)


Updated 21.08.2017