[Tobi]

Hey everybody,

I have the feeling that solving solid mechanics equations making me crazy and every day I feel a bit more stupid. What I realized and verified now is that if I have an equation and use a relaxation factor of 1, I get different calculation times than using no value (or a negative one):

• Using no relaxation will solve my validation case within 0.77s
• Using a relaxation factor of 1 will lead to an time increase of 32s

I checked the code and I know why this behavior take place but I have no idea why the code is implemented like that:

• In the fvMatrix.C source file we will find the code for relaxing the equation

1. We store the old diag matrix
2. We calc the sum of the off diagonals
3. We add the maximum contribution of the boundaries to the diag
4. We check if the sum of the off diags are larger than the diag itself, if so, put the off diag value to the diag
5. We relax the diag (increasing all values by dividing with relax factor)
6. We remove the minimum contribution of the boundaries to the diag
7. We put the relaxed stuff to the source

The solutions are the same but the time is different. I am asking why we add first the maximum and then remove the minimum of the patches to the diag. Does anybody have any references about that? In Moukalled et al. [21] we will find some explanation about the function and that it is the Panktar implicit relaxation but there is no information about the constraints based on the boundaries. Finally, I wanted to point out that using a relax factor of 1 changes my matrix based on the patch handling and I would like to know why we are doing it in that way or in general the information behind it I guessed that using 1 will not influence my matrix system as I could see from the equations that are given e.g. in [21]. I cannot believe that this is a bug. I mean first adding the maximum and then removing the minimum. Because if it would be twice the maximum, using the relax factor of 1 will lead to the solution I would expect (no influence). But I really sure that I am wrong here and there is a meaning about using first the maximum and then the minimum.

Hope some expert in that topic can tell me whats going on here.

else
{
    // For non-coupled boundaries add the maximum magnitude diagonal
    // contribution to ensure stability
    forAll(pa, face)
    {
        D[pa[face]] += cmptMax(cmptMag(iCoeffs[face]));
    }
}
.
.
.
after relaxing
.
.
.
else
{
    forAll(pa, face)
    {
        D[pa[face]] -= cmptMin(iCoeffs[face]);
    }
}

If someone would agree that it should be cmptMax in the second part too, I will make a bug-report.

[arjun]

I am of opinion that using no relaxation and relaxation factor of 1 shall produce same results.
So i personally would classify it as a bug.

Second about implicit relaxation, Milovan's (Peric) book I think mentions the derivation.

I bring your attention to the fact that this derivation is based on the assumption that one solves the linear system fully ie after the relaxation whatever system you get you are supposed to solve till machine precision then you get that final form of formula that is being used. This part is important because often they are not solved to machine precision. In some cases it does affect how yyou get solution but it is rarely mentioned in any text.

[Tobi]

Hi Arjun,

thank you for your reply. I have to mention (based on the last sentence you made - and I hope I got it correct) that I get the same solution but with much more effort. This is caused by the treatment of the already mentioned addition of the maximum and removing of the minimum components. For a scalar it should not matter (same values) but for a vector quantity it matters. The problem in the displacement calculation is the special treatment of the boundaries which will in fact influence the matrix with relax factor of 1.

By the way, Moukallet et al. [21] also derived the equations and based on that, and other literature, the relaxation factor of 1 should give the same matrix (as you confirmed). I will make a bug report but I guess Henry will just close it and mention that it is not a bug Thank you Arjun.

[Tobi]

I have no Idea but as far as I got it, Henry confirmed that the matrix does not has to be equal for vector matrices. Bug report closed. Relax 1 for vector matrices really influence the matrix itself (not the result).

https://bugs.openfoam.org/view.php?id=2522

[alexeym]

Hi,

Different times are the result of these lines (https://cpp.openfoam.org/v4/a00885.h...f646011c212e):

if (alpha <= 0)
{
    return;
}

If you do not set relaxation factor alpha is equal to 0. So in both cases (no relaxation information or negative value) relaxation code is just skipped. While if alpha == 1, relaxation code is executed, though with no visible results (except for execution time), since D0 == D. I think, in terms of the Foundation it is called feature not a bug.

[Tobi]

Hi,

I agree with alpha <= zero that we skip it. But alpha = 1 will change a vector matrix if the components at the boundaries differs. I would agree that alpha = 1 will not affect the matrix, like you said and arjun however this is wrong for vector matrices and I will show why. At the moment there are two cases in which we will modify the matrix while alpha = 1:

a) The matrix is not diagonal dominant, therefore the diagonal element which is not dominant will be replaced by the sum of the off diagonals:

  621     forAll(D, celli)  
  622     {
  623             D[celli] = max(mag(D[celli]), sumOff[celli]);
  624     }

This will hold for all types of matrices.

b) The matrix represents a non-scalar matrix and the boundary values in the direction differs. Why? Because first we add the magnitue of the maximum component:

  558                 // For non-coupled boundaries add the maximum magnitude diagonal  559                 // contribution to ensure stability
  560                 forAll(pa, face)
  561                 {
  562                          D[pa[face]] += cmptMax(cmptMag(iCoeffs[face]));
  563                 }

and then we subtract the minimum component:

  648                 forAll(pa, face)
  649                 {
  650                         D[pa[face]] -= cmptMin(iCoeffs[face]);
  651                 }

If we would add and subtract the same values, alpha = 1 will not affect the matrix in case b) but in case a). But a) is wanted to improve the matrix.

So case a) is a good thing but in case b) I am not sure if this is a mistake or there is a meaning behind it. That was my question and I am sorry if I did not make that statement clear enough. Thus, alpha = 1 will not only affect the time.

[arjun]

I also knew he would mark it as feature and not bug because he intended that behaviour. However alpha = 1 shall be equivalent to no relaxation.

Assuming this is not a bug as he said, i would open another thread asking for explanation why alpha = 1 and no relaxation performance results are different for inconsistent behaviour. I know it would bug him but the behaviour shall be consistent with no relaxation which is not.

Anyway, when matrix is not diagonal some softwares make it diagonally dominant. I am of opinion that actual matrix shall not be touched in a way that is not consistent with literature. Anything special behaviour user shall decide to switch on or off.

Some equations do not like diagonals be touched in anyway as this is the case.

[lengweee]

Hi, Tobi,

My comments on openfoam handling relaxation is as follows.

First,
foam stores the linear system in the way that internal and boundary info are separated, e.g., instead of store A, f to solve Au = f,
it stores A_i, f_i for internal info,
and A_b, f_b for boundary info,
and solves (A_i + A_b) u = f_i + f_b.
See

https://www.cfd-online.com/Forums/openfoam/167028-extracting-matrix-data-fvvectormatrix.html#post703887

for export the linear system.

The reason the matrix is stored separately is to save storage and computation,
e.g. the internal matrix is the same for the three component of velocity of incompressible flow.


Second, generally, the relaxation applying to A * u = f is to modify D(A) to be

D(A)_new = D(A)/alpha

, or

D(A)_new = max{D(A), sumOff} / alpha,

, where D() is diagonal, and update RHS as

f = f + (D(A)_new - D(A)) * uold

. The maximum is taken to make diagonal dominace.

However, for the separated case, foam choose other ways to do it.

First way, from older version, change both internal and boundary info.
1. on boundary,
update diag with

    D(A_b)_new = | D(A_b) | / alpha,

update rhs with

    f_b,new = f_b + ( D(A_b)_new - D(A_b) ) * x

2. in interior
update diag with

    D(A_i)_new = max{ D(A_i), sumOff - D(A_b)_new }

or in another way

   D(A_i)_new + D(A_b)_new = max{ D(A_i) + D(A_b)_new, sumOff }

which is an approximation of

    D(A)_new = max{D(A), sumOff},

update rhs with

    f_i,new = f_i + ( D(A_i)_new - D(A_i) ) * x

Second way, from newer version, change only internal info.
Add positive boundary contribution to diag

     D = D(A_i) + |D(A_b)|    ... (*)

relax diag

     D_new = max{D, sumOff} / alpha

update interal diag

     D(A_i)_new = D_new - D(A_b)  ...(**)

update rhs with

     f_i,new = f_i + ( D(A_i)_new - D(A_i) ) * x

To achieve diagonal dominance,
in equ(*), max is applied for components,
in equ(**), min is applied for components.

Comment #1, if all components has same boundary type, the A_b is the same for all components.

Comment #2, ONLY when A_b is POSITIVE that the second way is equivalent to apply relaxation to A and f directly.

Comment #3, ONLY when A_b is POSITIVE and A is diagonal dominance that the second way is compatible to original problem. Otherwise the relax scheme is solving ANOTHER problem !

Comment #4, way 2 seems more concise.