[chriss85]

Here are my results using the settings from the case I posted. Unfortunately the image gets compressed quite much here, but the results should be visible.

Do you think that higher precision can be reached somehow?
The results seem to converge for a relaxation factor of 1. This may or may not mean that these results are the most correct. Unfortunately in a real case I need relaxation to archieve convergence. I should probably also test if this can be avoided by reducing the time step size. I need some kind of estimate on the convergence behaviour of the iteration as a input signal for finding the proper time step size though.

[Tobi]

Hi Chriss,

you are welcomed for discussion. Also I learn always new things and repeating is very good for keeping the stuff in mind. Okay as far as you told, you get the point of relaxation, final etc. but in the time I made this PIMPLE blog, I outputted the relaxation factors and indeed the relax factors that were using in the last final iteration kept 0.x if you do not set the Final entry (well I am still not 100% sure if I mix up sth. because you showed the two lines with the default values). I will check this again to be sure to give you a correct answer.

I agree, that if you are using no residual control, you should have the a high accurate solution. The best would be to have no relaxation and using a very small time step to be stable.

To your relaxation factors and zero. I do not agree with you. My college also told me that the equation is:


So if you set alpha to zero, you should not go on in the solution.
I check this!

PS: what is the y-axis in your plot?

PPS:

    if (alpha <= 0)
    {
        return;
    }

As you see, if you insert some non-physical value for alpha it is not running the function. The function is getting quited and no relaxation is used. That means, that it is the same like no relaxation.
I will output the relaxation factors and give a feedback.

Also you can see it a bit below:

// ... then relax                                    D /= alpha;

Alpha = 0 would be give an floating point exception. And alpha = 1 does not affect the matrix D

19.06.2019... the statement above is wrong. Alpha = 1 affects the matrix. It makes it at least diagonal dominant.

[chriss85]

We need to differentiate between field and equation relaxation. Field relaxation is using the formula you mentioned and keeps the old value for alpha=0 and only uses the new value for alpha=1 (so no relaxation in this case). The code I linked to is about equation relaxation though (this is why it's in fvMatrix.C). If you look a bit before D/=alpha you will see that it's already modifying the matrix to make it diagonally dominant, even when alpha=1, so alpha=1 is not necessarily the same as not relaxing. 0 on the other hand is because of the return statement.

The y-Axis values are in SI units, so velocity in m/s, pressure in Pa (N/m²), Temperature in K.

Maybe we should check the final relaxation factors again, just to make sure. Something might have changed since the time when you tested this or maybe there
s something else in the code I didn't consider.

[Tobi]

Hi all, dear Chriss,

I checked the relaxation function and here is the result:

### Case 1 ###
Using the following relax dictionary:

relaxationFactors
{
    fields
    {
    }
    equations
    {
        e       0.3;
        U       0.6;
    }
}

will relax all your iterations but not the final iteration will not be relaxed due to the fact that there is no setting for UFinal/eFinal. This can be checked here:

https://github.com/OpenFOAM/OpenFOAM...vMatrix.C#L685

The output of the last two pimple iterations are:

PIMPLE: iteration 8
UEqn.relax() start
In function relax() --> name:U
In function relax() --> psi_.mesh().relaxEqn(): 1

    From function fvMatrix<Type>::relax(const scalar alpha)
    in file /home/shorty/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/fvMatrix.C at line 534
    Relaxing U by 0.6
UEqn.relax() end
DILUPBiCG:  Solving for Ux, Initial residual = 5.103798e-06, Final residual = 1.1441873e-09, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 3.6895691e-06, Final residual = 9.9612322e-10, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 5.1905722e-06, Final residual = 8.8635583e-10, No Iterations 1
In function relax() --> name:e
In function relax() --> psi_.mesh().relaxEqn(): 1

    From function fvMatrix<Type>::relax(const scalar alpha)
    in file /home/shorty/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/fvMatrix.C at line 534
    Relaxing e by 0.3
DILUPBiCG:  Solving for e, Initial residual = 2.0190539e-05, Final residual = 8.822835e-10, No Iterations 1
DILUPBiCG:  Solving for p, Initial residual = 0.00029119078, Final residual = 5.8967572e-07, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.6310611e-06, global = -1.3652136e-06, cumulative = 2.4972179e-06
DILUPBiCG:  Solving for p, Initial residual = 6.0992289e-07, Final residual = 5.3852368e-09, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.4900303e-08, global = -1.2560377e-08, cumulative = 2.4846575e-06
PIMPLE: iteration 9
UEqn.relax() start
In function relax() --> name:UFinal
In function relax() --> psi_.mesh().relaxEqn(): 0                 // - no relaxation
UEqn.relax() end
DILUPBiCG:  Solving for Ux, Initial residual = 4.0998722e-06, Final residual = 4.0930334e-12, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 3.0195606e-06, Final residual = 2.0841141e-12, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 4.2793435e-06, Final residual = 3.3578357e-12, No Iterations 2
In function relax() --> name:eFinal
In function relax() --> psi_.mesh().relaxEqn(): 0                 // - no relaxation
DILUPBiCG:  Solving for e, Initial residual = 4.8759542e-05, Final residual = 2.5880915e-11, No Iterations 2
DILUPBiCG:  Solving for p, Initial residual = 0.00070397614, Final residual = 5.3164087e-11, No Iterations 6
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.4674538e-10, global = -1.4128755e-10, cumulative = 2.4845162e-06
DILUPBiCG:  Solving for p, Initial residual = 1.7891292e-07, Final residual = 4.966156e-12, No Iterations 3
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.3721562e-11, global = -9.4849471e-13, cumulative = 2.4845152e-06
PIMPLE: converged in 9 iterations
ExecutionTime = 0.51 s  ClockTime = 0 s

End

### Case 2 ###
Using the following relax dictionary (typical error):

relaxationFactors
{
    fields
    {
    }
    equations
    {
        "e.*"       0.3;
        "U.*"       0.6;
    }
}

The output of the last two pimple iterations are:

PIMPLE: iteration 8
UEqn.relax() start
In function relax() --> name:U
In function relax() --> psi_.mesh().relaxEqn(): 1

    From function fvMatrix<Type>::relax(const scalar alpha)
    in file /home/shorty/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/fvMatrix.C at line 534
    Relaxing U by 0.6
UEqn.relax() end
.
.
.
In function relax() --> name:e
In function relax() --> psi_.mesh().relaxEqn(): 1

    From function fvMatrix<Type>::relax(const scalar alpha)
    in file /home/shorty/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/fvMatrix.C at line 534
    Relaxing e by 0.3
.
.
.
PIMPLE: iteration 9
UEqn.relax() start
In function relax() --> name:UFinal
In function relax() --> psi_.mesh().relaxEqn(): 1

    From function fvMatrix<Type>::relax(const scalar alpha)
    in file /home/shorty/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/fvMatrix.C at line 534
    Relaxing U by 0.6
UEqn.relax() end
.
.
.
In function relax() --> name:eFinal
In function relax() --> psi_.mesh().relaxEqn(): 1

    From function fvMatrix<Type>::relax(const scalar alpha)
    in file /home/shorty/OpenFOAM/OpenFOAM-2.3.1/src/finiteVolume/lnInclude/fvMatrix.C at line 534
    Relaxing e by 0.3
.
.
.
PIMPLE: converged in 9 iterations
ExecutionTime = 0.48 s  ClockTime = 1 s

End

You can see, that the last iteration will still use the under relaxation factors. And that will not be time conservative.

And in Case three you can imagine what happen if you introduce the Final relaxation factors and set all to 1. It should be the same as the first one but the second one is wrong.


As you mentioned before, that we have to make a difference between the field and equation relaxation function, I agree.
In equation we are using matrixes and the diagonal elements. D0 stores the unrelaxed diagonal and D will be relaxed. At least D will be relaxed with D /= alpha and at least the sources will be calculated again:
S += (D- D0)*psi_.internalField();

If Alpha = 1 --> S keeps the same.
I did not go through the whole diag matrix calculation but it should be.
I will check the last statement and give a feedback.


PS: SI for Temperature in K --> 50.000 K is that correct?

[chriss85]

Yes you're right with this, wildcards/regular expressions should not be used in the relaxation dictionaries.
By the way, you can use the debug mechanism for outputting things like this if they're already in the code.

See here: https://openfoamwiki.net/index.php/H..._from_OpenFOAM

and enable the debug flags for fvScalarMatrix / fvVectorMatrix.
You won't have to recompile anything then.

So, back to the topic, any ideas about the differences I get with different relaxation factors?

[chriss85]

I've also been thinking if it's maybe more reasonable to use smaller time steps and no relaxation at all. If the problem becomes stable under these conditions it might be the better option. The only question is how this would influence the calculation time. One would probably save some PIMPLE loops per time step but have more time steps in total which probably more than compensates the savings.

I could write a control loop that modifies the time step size according to a target number of PIMPLE iterations while keepinga CFL number as an upper limit for the time step size. Do you think this could be a sensible approach?

[Tobi]

Originally Posted by chriss85

By the way, you can use the debug mechanism for outputting things like this if they're already in the code.

See here: https://openfoamwiki.net/index.php/H..._from_OpenFOAM

Thanks, I already know that, but if you want additional output you have to hard code it
Back to the topic,...

I am not sure if in special cases D0 and D differs if you have alpha = 1. I think not.
In your case it`s the same. I ran the case and outputted the matrix elements for D0 and D in the final iteration of one pimple loop:

~$ diff oldMatrix newMatrix
~$

Both matrixes are identically, see also in the plot.
Anyway I am not so good in mathematics than my college and unfortunately he is in holiday to ask him about that.


Back to the origin topic:
I checked out your case with different settings.

1. no relaxation
2. no relaxation | other dt like in 1.
3. relaxation 0.6 | same dt like in 1.
4. relaxation 0.2 | same dt like in 1.
5. relaxation 0.2 | same dt like in 2.

I solved until I reached 1e-8. The problem is a bit complex and for me not so easy to handle because I have never used sonicFoam and such cases (a lot of things happen in between a short dt) where sonic or supersonic flow occur.

Well lets try to get some good information out of it.
Let`s compare simulation 1, 2, 4 and 5 after 1e-9s. For comparison I made the same pressure contour interval for each case (its figure p1.png).

• Simulation 1 vs 2 is clear for me. In 2 we have a bigger time step and therefore we lose information. The extreme thing in that case (I think) is that at the very first micro seconds a lot of things happen, which needs a good time resolution to get the information. Otherwise, it is lost.
  
  In this case #2 ran one iteration and case #1 ran 20 iterations.

• Simulation 2 and 3 are looking similar but I have no idea why. Normally 3 should be more like 1 due to better time accuracy and 4 should be like 2.

If we compare the solution after 8e-9s (p2.png) we can see that the solution is really equal to each other. I think that depend mostly due to the fact that in the early stage you need a very good time accuracy to get the direction/path of the information.




In figure p3 and p4 you can see the same for t=1e-9s and t=8e-9s but this time with the three relaxation. I have no time to investigate more into that topic. A bit strange but maybe my abort criterion is too low (1e-8) or I made some very simple mistake.

Hope it will help you in some way.

[Tobi]

Originally Posted by chriss85

I've also been thinking if it's maybe more reasonable to use smaller time steps and no relaxation at all. If the problem becomes stable under these conditions it might be the better option. The only question is how this would influence the calculation time. One would probably save some PIMPLE loops per time step but have more time steps in total which probably more than compensates the savings.

I could write a control loop that modifies the time step size according to a target number of PIMPLE iterations while keepinga CFL number as an upper limit for the time step size. Do you think this could be a sensible approach?

I think it's a good start. As I mentioned, the time resolution is necessary to get the information path (I only read about it, because I never used sonic/supersonic flows). Normally the Pimple loop should use the Co criterion but as far as I know, sonicFoam does not use it.

At the moment I can not help you anymore in that topic. Sorry for that.

But keep in mind that we are always using some models and interpolation schemes etc.
It's funny because if you replace a turbulence models, grid, schemes etc. you will always get other results.
That's why its important to understand the mathematics behind all the stuff (but I am still a Ph.D. student and learning all the day )

Kind regards
Tobi

[chriss85]

Adding adaptive time step sizes based on Courant number to sonicFoam is quite easy and one of the first things I did
However, sometimes this adaption is not enough for me and I need some more adaption based on other indicators. These can be field convergence behaviour (is it oscillating during the iteration?), number of required PIMPLE loops for convergence, residual behavior, ... . The problem lies in finding a good indicator and a good control algorithm. Of course we can also resort to using smaller Courant numbers, I just feel I'm wasting some performance then.

A further complication is that I have figured out that it is beneficial to calculate the radiation only every few iterations, because it is quite costly and the loop converges faster (calculation time) when it isn't calculated every step. This can make differences between two iterations unreliable because the convergence behavior might show steps.

[Tobi]

Originally Posted by chriss85

Adding adaptive time step sizes based on Courant number to sonicFoam is quite easy and one of the first things I did

In the previous post I wanted to say that you have to implement it yourself
Everything is still included so I know that its not a hard task

So good luck and keep the community up to date.

[chriss85]

I have performed a calculation using my sonicFoam-based solver with deltaT=x and deltaT=0.5x with x being some arbitrary value. I run the iterations until convergence is archieved (i.e. residuals don't improve much more).
Runtime roughly doubles and number of PIMPLE iterations stays more or less constant. This means that the number of required iterations doesn't correlate much with the time step size.

I also tried to run it using a combination of absolute residuals and relative residual values. In this case the runtime was a bit lower than the double time, but still significantly higher.

So all in all I think it's not worth it to use lower deltaT values to save PIMPLE iterations. It may be needed for stability though to avoid the need for relaxation.

[danny123]

Hello Chriss,

Since you asked for the number of Pimple loops (I forgot the thread, but this one sounds good too):

Number of pimple loops for residual control is ideally around 2 to 3. Beyond this value, I get divergence of the p field in particular. This value allows stable implicit deltaT control for a given initial residual of p_rgh of 5e-04 (this is a number using my own method, not the OpenFoam residual number). I think increasing that initial residual number does not lead to a different required PIMPLE number, but the required deltaT is different.

So, the method works. I get a tremendeously low residual number for alpha1Eq, p_rghEq is more difficult. Still the required deltaT is very small (10^-7). I am now looking into the inner loops (PISO or NonOrthogonal Corrections) to improove that situation. Unfortunatly it picks up the wrong loop doing so. This is why I have opened a new thread on that topic.

Regards,

Daniel

[chriss85]

I wasn't able to do something similar in sonicFoam unfortunately. When I use half the value of deltaT, the required number of PIMPLE iterations stays nearly constant for me

Also, I would be more than happy if I could use timesteps of dt=1e-7 in my case
I'm usually two orders of magnitude lower, but the total simulation time is also on the order of 1e-5 to 1e-4.

[Raza Javed]

Originally Posted by Tobi

Hello Chriss,

some hints to your case and underrelaxation:

• you have to specify two relaxation factors for one variable (e.g.)

 
U    0.4
UFinal 1

• if you don't use UFinal you will always have 0.4 and then the final iteration in the pimple loop is also done with 0.4 (as I remember correctly). That means that you have to run xxxx outer iterations to be sure your solution is converged (but that is not the sence of underrelaxation)

• using PIMPLE you should define residialControl parameter, otherwise you will run 200 outer iteration (in your example). That is very senceless and expensive.

• Additionally you should always define the solution directory for relaxed and final iteration. Here you should define relTol 0.1 for the solver and in the final directory relTol 0

You can download the case here:

www.holzmann-cfd.de/cfd-online/cases/cfdOnline_SonicFoam_chriss.tar.gz

I did not check the results till its converged but the results should be the same. Please notice that with underrelaxation you will reach a bigger timestep and therefore, you can loose timer information.

Kind regards
Tobi

Hi,


I have one question related to relaxation factor.


I am using chtMultiRegionSimpleFoam, and my openfoam version is 4.1.


My geometry has one box region, and inside that box region, there is one fluid region and three heaters (making heat source using fvOptions)



My fluid flow is Turbulent, and I am using kEpsilon Turbulence model.


I have compressible fluid flow.


Now, I have set all the boundary conditions for p,U,T,k,Epsilon,nut,alphat.


Now, the physical expectation of the case is that, I am putting 300K temperature at the inlet of the fluid, and then I want to see different temperature at the outlet of the fluid region. (temperature change due to heaters)


Now, I am seeing this effect.



When I RUN the simulation for 50 seconds, then from the temperature difference between intlet and outlet, and the heat capacity of fluid(water) I can find that how much power this fluid is throwing out of the system.


And after hand calculation, i checked that this is the same power that I am putting on three heaters using fvOptions.


it means that the simulation results and the hand calculation results are same.



Here comes the problem:


When I RUN the simulation for longer time (say 100 seconds), the temperature at the outlet starts increasing, and this makes the greater temperature difference between inlet and outlet. and then the hand calculation goes wrong.


If I increase the simulation time further, the temperature at the outlet rises even more.


Now, I don't know that what could be the problem here?


Do I need to change the relaxation factor to make the system stable? OR something else?


I shall be very thankful if you can help me in this. I would be happy if you need any further information about my case.


Thank you

[Tobi]

I guess it is related to the boundary condition, not the relaxation factors but out of the box I cannot give you an correct statement.

[Raza Javed]

Originally Posted by Tobi

I guess it is related to the boundary condition, not the relaxation factors but out of the box I cannot give you an correct statement.

Thank you so much for your reply.


Yes, it was the error of boundary conditions.


I have one question that is related to Paraview. I know this thread is about relaxation factor, I asked this question in post-processing part, but couldn't get response, so I am asking you here.


In my geometry, I have one fluid region with inlet and outlet. I want to plot the graph of temperature and time at the outlet, that how the temperature is changing with respect to time at the outlet.


I tried using extract block filter, with that I can select the outlet but I am unable to obtain the temperature vs time plot at the outlet. can you please help me in this?


And if I cut my fluid region in half, and want to see the same graph on the surface I got after cutting, how can I do that?


Thank you so much for your help.