[chriss85]

The system I'm solving has a very wide range of dynamics. It's based on sonicFoam, so this means the PIMPLE scheme is employed.

Relaxation is used to stabilize the system by mixing the current result with the previous result (or doing something similar to the matrix).
So generally, low relaxation factors mean stable but slow convergence, high relaxation factors are potentially instable but might converge faster.
Best convergence is reached if the relaxation factors are chosen in such a way that no (or only small and decreasing ) oscillations occur in the resulting field.

I am wondering if it would be a good idea to adaptively adjust the relaxation factors based on field probes behaviour. For example, one could watch the minimum/maximum field values and see if they oscillate during the convergence loop. If they do, some relaxation factor is probably too high.

Now, the oscillating field isn't necessarily only coupled to the relaxation factor of the equation for this field. Because of the coupling between different equations the relaxation factor for one equation may influence the oscillatory behaviour of another field.

If one would adapt all relaxation factors at once, it would be possible to get rid of oscillatory behaviour, but at the price of sacrifizing speed because some relaxation factors are not correct.
Maybe it is possible to adapt each factor separately, but only if a relaxation factor causes oscillations in its own field and possibly others, but not only in other fields.

Any better suggestions here?

Another problematic that comes up here is finding the correct breaking conditions for the convergence loop. I'm currently breaking when the initial residuals haven't decreased by a certain value over a certain number of steps. I feel this is more reliable rather than setting fixed absolute values or relative values to the first loop iteration residual. However, if the relaxation factor changes, the speed of convergence will also change, and the breaking condition may either get too fast or too slow. Any ideas on this issue?

Generally, I believe that with properly adapted relaxation factors and breaking conditions a speedup in a typical case of about 2x should be realistic, so this should be a problem worth pursuing.

I welcome any ideas and comments to my thoughts

[danny123]

Since you started a new thread, here my experience, which is not huge and definetly not very positive.

I do relaxation of the PCG solver (and other matrix solver) in order to correct if the solution diverges. This can save calculation time sometimes as long as you do not overdo it. If the relaxation factor is small, the whole PCG scheme is messed up and calculation time increases exponentiallly.

Since my problem is transient, relaxation in time is not permitted. It means that the time to get to convergence is longer. But this is unphysical since this is not time of the correct Navier-Stokes equation anymore.

One can relax in between, meaning you get a solution closer to the final one, then switch to relaxation 1. This is what OpenFoam does by default, but this is unstable in my case.

I think the whole thing is related to whether the new pressure field influences a lot the flow compared to old onces. This is the case for my problem. This means relaxation cannot be used, also not the moment predictor. As has to be switched off.

OpenFoam uses a relative error measurement for the matrix solver (even if it is called abstol). The error is put in relationship to the error when the solver starts. The error at the start depends on how close you are to the solution at beginning of your iteration. If this error decreases quickly this means your convergence rate is quick. One assumes that for low changing convergence rates your are close to the final solution. But this you cannot know, your solver might just be slow.

Regards,

Daniel

[chriss85]

Are you saying that you can not use relaxation for your problem?

I actually need to use some relaxation for stabilizing my convergence in a time step, however the amount that is needed can change during a simulation, which is why I was talking about controlling the relaxation factors automatically somehow.

Do you (or anyone else) happen to know why sonicFoam doesn't relax the pressure, like other compressible solvers do? I get some fluctuations on the pressure fields which I would like to relax, but I don't know if there is a specific reason this is not done?

[chriss85]

I just noticed a huge problem in my solver. When I use different relaxation factors, my results are affected, even though the solution appears to be converged (by looking at residuals and maximum field values). Any idea what could be causing this? Attached are some plots that demonstrate this problem with relaxation factors for energy and velocity of 0.1 and 0.3.

[danny123]

You have to set the final relaxation to 1 (Ufinal or pfinal) to ensure that at least at the end of the iteration, you comply to the Navier Stokes eq.

To ensure consistency of an iterative equation, it is necessary to relax the change within a guess value. Let's say you have a Fix point iteration. If you relax such a method, the result is wrong. You can speed up the iteration or slow it down by choosing a relaxation factor other than 1. But at a certain point (usually when the convergence rate becomes flat), you need to switch back to 1. This is also how OpenFoam relaxation works.

The alternative is to relax the difference between the previous and the current estimation at relaxation 1. This is how e.g. the sekant procedure is relaxed. Basically you relax a function, which is zero when you have found the solution instead of the iteration value.

Another point may be that the so called normfactor is different between iterations that start with a different relaxation factors. This normfactor ensures that the residual is 1 at the beginning of the iteration. If you apply a different relaxation factor, I am not sure that you get the same normfactor, which means I am not sure that you even can compare residuals for different relaxation factors. This is easy to test by plotting that normfactor out (it is within the code of the matrix solver e.g. PCG.C).

Regards,

Daniel

[chriss85]

I have just verified that it is not required to specify the UFinal, eFinal, .. settings manually, as the code uses default values (1 for explicit field relaxation, 0 for matrix relaxation). Interesting fact: relaxing a matrix with a relaxation factor of 1 is not equal to doing nothing. Instead, it makes the matrix diagonally dominant. I believe this can help the convergence of the matrix solvers sometimes if the matrix is ill conditioned. A value of 0 is handled specifically in the relaxation to do nothing.

Also, when we talk about explicit field relaxation I don't believe that you will get a correct result in any case, since the whole algorithm is an iterative procedure. Where you arive with the result depends on the previous iteration field which is used for relaxation. This is determined by calling field.storePrevIter(). If you do this directly before solving the corresponding equation it should be fine. Even if you relax with a lower factor than 1 your result can be (within the approximative limits of the convergence loop) correct, if the loop has converged good enough already. On the other hand, if your system were to be extremely unstable, using a relaxationFactor of 1 might blow it up in the final step. I don't think this is a problem in real applications though.

I'm also not quite sure what you mean with the two different relaxation methods you mention. In my understanding, field relaxation does U_relaxed = U_(n-1) + f*(U_n - U_(n-1)), where U_n is the current field value and U_(n-1) is the value stored by U.storePrevIter(). This is basically equal to relaxing the difference between the two fields.

I will try to look into the normfactor, as I haven't looked much into the matrix solving code itself yet. However, the problem I mentioned in the previous post is also very visible in the maximum field values, which are completely different when the relaxation factors are changed. It's basically the same calculation but only with different relaxation factors. It appears to me that this means that the iteration doesn't have a single fix point and the converged value is not unique???

[danny123]

Well, this blow up on the final application did happen in my case. So it is ideed real at least for my application. I think it depends on how strong the coupling between velocity field and pressure is. This is also why the momentum predictor does not do any good in my case.

I do not use the relaxation of the matrix that is implemented in OpenFoam, but I built up my own code, which works differently.

In my case, the pressure field and U field are periodically changing (due to the agitation of the stirrer). When I used different relaxation factors in the past, the wave form were similar, but average values differed. Everythings works fine if I have a continuous steady state problem. This is why I gave up on this.

Relaxing is a very brute method manipulating iterative solvers. They are typically designed as a linear extrapolation of the previous guess value using a kind of Newtonian method. For some problems overshooting means that your solver does not find back since the function makes some strange turn. In this case you restrict that overshooting by relaxing.

Regards,

Daniel

[chriss85]

If you see different results with different relaxation factors, I would expect that:

a) Your solution wasn't converged
b) You have some kind of code that doesn't converge to a unique solution

I have now tested the case I was simulating with vanilla sonicFoam and I also see a dependence on the relaxation factors. This is really troubling
I'm going to try to simplify it as much as possible until I can identify the problem.

[chriss85]

I have uploaded a test case for sonicFoam where I get different results with different relaxation factors. Could you run this case with relaxation factors for e and U of 0.1 and 0.3 and check for differences? I see different maximum field values even though the solution appears to be converged.

http://s000.tinyupload.com/index.php...96535611154888

Do you think these low relaxation factors lead to deviations in the equations being solved? I tried a tutorial case with the same relaxation factors where I got the same result no matter what relaxation factor was used.

[danny123]

Maybe you take a look at this:

http://holzmann-cfd.de/cfd-online/OpenFOAM.pdf

As a number of other literature leads, it states that relaxation is not time conversative. Look at case C.

I have not verified these statements, neither have I the time to do so. I think relaxing is done to avoid oszilatory behavior of steady state problems. You use the new solution and correct it with the previous one using some factor. So, if your solution oszilates, the wave amplitude is flatened, this is what this procedure was built for.

If you have a transient problem, you just defer the correct solution calculation backwards. The solver has a reason to diverge, which is not fixed. The resiuals flatens out earlier, but this does not mean necessarily that you are any closer to a final, best solution.

Regards,

Daniel

[chriss85]

I agree with you that relaxing over different time steps will lead to a solution which is not time conservative. That is why relaxation is being done inside the iterations of a single time step. If you use the previous iteration result as a basis for the field relaxation it should simply stabilize the convergence (and possibly slow it down). If the code looks like this then it shouldn't be a problem, as long as you run the iteration in the time step until convergence is reached:

while(runTime.loop())
{
  while(!converged)
  {
    fieldA.storePrevIter();
    fieldAEqn.solve();
    fieldA.relax();

    fieldB.storePrevIter();
    fieldBEqn.solve();
    fieldB.relax();
  }
}

If you don't use enough iterations, the solution contains some contributions from the previous time step (from the first call to relax which uses previous time step field) and from the unconverged, wrong values. However when the system is stable enough it should be able to reach a converged solution. As long as there is only one unique convergence point this should then be independent from the relaxation factors. Only the number of required iterations shoud change.

If your code looks like this on the other hand, than the last iteration must not use field relaxation:

while(runTime.loop())
{
  fieldA.storePrevIter();
  fieldB.storePrevIter();
  while(!converged)
  {
    fieldAEqn.solve();
    fieldA.relax();

    fieldBEqn.solve();
    fieldB.relax();
  }
}

There are possibly other factors that lead to wrong results, such as too low solver tolerances, bad mesh, possibly boundary conditions.
I have tried lowering the solver tolerances by a few orders of magnitude and the mesh is ok according to checkMesh. I have only used adiabatic walls with no slip conditions, so nothing special. Any other ideas?

[danny123]

This looks ok to me except that the loops are imbedded not consequetive. But this is a detail. Did you look at conflicting iteration loops?

In your case, I looked at your solver quickly (I am by no means expert on this one). Maybe you look at those areas:

- your thermodynamic coupling may oscilate or rhoEq.
- the moment predictor may be in conflict with the reconstructed field after pEq
- the phi field calculation may diverge.
- your energy balance may introduce oscillation (not very typical though)

Regards,

Daniel

[chriss85]

I have done some more tests, namely by initializing the mesh to different field values. When I use much lower temperatures and pressures the results appear to be better, so I suspect some numerical problems somewhere. I'm going to try some more things like setting temperature calculation tolerance and iteration number higher. If this doesn't help than I fear that this is a problem/limitation of the sonicFoam solver and I will have to find some other code which is better suited to high speed, temperature and pressure flows.

[Tobi]

Dear Chriss,

I checked your case but I dont get the point of your simulation because the file 0 is missing. Also you don't set any pressure or velocity BC ? For U everything is (0 0 0)... Actually I dont know what your case is used for and if its a closed geometry, some geometry with outlet/inlet etc.

If I simulate you case for one iteration it seems that your BC are not really describing what you wanna do.

For example:

• no pressure
• no velocity
• no fvOptions

Due to your BC I see that you named one pach "outlet" and you have some ignition patch. Can you describe the case and add the 0 folder (I already created one out of your latest time step but there are no special BC).




First be sure to have a mathematical correct set of equation that describe your problem.

[chriss85]

I'm really sorry about that, I must have missed something when packing up the case. Here is a fixed version I just tested, just execute sh AllRun: http://s000.tinyupload.com/index.php...34340727958426

The boundary names are leftovers from my real application, please ignore them. In this case I specify a higher pressure/temperature in a specific region and only use closed walls and adiabatic boundaries to ensure there is no interference from the boundaries. I also still have a velocity field from the real case, but it can be anything really.

I'm going to post results for this case very soon.

[Tobi]

Hi Chriss,

you made a lot of mistakes using underrelaxation. So I think you did not get the point in my blog

I will run your simulation and give a feedback.

[chriss85]

I don't think so. I specifically did not want to use residualControl here if that's what your meaning, to make sure that with every tested relaxation factor the solution is converged to the best possible value. It's clear that I'm sacrificing performance here, but this is intended to compare precision as good as possible. In my real case I use a much more elaborate residual control (I have posted this recently on the forum). The settings used here are only meant to simplify the problem as much as possible and use the best possible precision.

[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) - 19.06.2019, the final relax factor if not defined should be 1 as default not the same as 0.4.. 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

[Tobi]

Originally Posted by chriss85

to make sure that with every tested relaxation factor the solution is converged to the best possible value.

So you did not get the point
You are using relaxation wrong, as mentioned above.

[chriss85]

I perfectly do get the point and I use residualControl in real world applications of course. However, using residualControl requires that you adjust your target residuals to the case in use. Depending on the system being solved, current time step and so on, different target residuals are needed. Relative values can be used as well, but what I did here was meant as simplification to make sure I don't add up tiny errors by possibly cancelling the iteration too early.

Regarding your other point, specifying UFinal relaxation factors (and similar) is NOT needed, because default values for this are already used in the code (check https://github.com/OpenFOAM/OpenFOAM...olution.C#L104 ).
Another thing to keep in mind is that when relaxing equations, using a value of 1 doesn't equal to not relaxing, instead a value of 0 needs to be used here (check https://github.com/OpenFOAM/OpenFOAM...vMatrix.C#L529 ). Relaxing an equation with a value of 1 ensures that the equation is diagonally dominant. I believe that in many cases this is equal to not relaxing, but it is not guaranteed.

Same thing goes for the solver settings in fvSolution here: Of course you are right that it is beneficial to use final settings for the equation solvers. Here again I have used very low tolerances everywhere to enforce best possible solution. This case is not meant to produce best performance but best precision. In fact I don't care about performance at all in this case, as it is only meant as a demonstrator for the influence of the relaxation factor.

By the way, thank you very much for your series of articles. It helped me very much to grasp the concepts when I started here