CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Information about equation relaxation (patch contribution)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By arjun
  • 1 Post By arjun

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 6, 2017, 05:55
Default Information about equation relaxation (patch contribution)
  #1
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,267
Blog Entries: 6
Rep Power: 39
Tobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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.

Code:
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.
__________________
Keep foaming,
Tobias Holzmann

Last edited by Tobi; April 6, 2017 at 08:09. Reason: small mistake - we remove and add the boundary stuff to the diag
Tobi is offline   Reply With Quote

Old   April 6, 2017, 07:44
Default
  #2
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 824
Rep Power: 22
arjun will become famous soon enough
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 likes this.
arjun is offline   Reply With Quote

Old   April 6, 2017, 08:17
Default
  #3
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,267
Blog Entries: 6
Rep Power: 39
Tobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   April 6, 2017, 12:03
Default
  #4
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,267
Blog Entries: 6
Rep Power: 39
Tobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   April 6, 2017, 17:06
Default
  #5
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,887
Rep Power: 34
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Hi,

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

Code:
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.
alexeym is offline   Reply With Quote

Old   April 6, 2017, 17:27
Default
  #6
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,267
Blog Entries: 6
Rep Power: 39
Tobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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:
Code:
  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:
Code:
  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:
Code:
  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.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   April 6, 2017, 20:05
Default
  #7
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 824
Rep Power: 22
arjun will become famous soon enough
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.
Tobi likes this.
arjun is offline   Reply With Quote

Old   August 25, 2018, 05:28
Default
  #8
New Member
 
wei leng
Join Date: Mar 2014
Posts: 3
Rep Power: 7
lengweee is on a distinguished road
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
HTML Code:
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
Code:
D(A)_new = D(A)/alpha
, or
Code:
D(A)_new = max{D(A), sumOff} / alpha,
, where D() is diagonal, and update RHS as
Code:
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
Code:
    D(A_b)_new = | D(A_b) | / alpha,
update rhs with
Code:
    f_b,new = f_b + ( D(A_b)_new - D(A_b) ) * x
2. in interior
update diag with
Code:
    D(A_i)_new = max{ D(A_i), sumOff - D(A_b)_new }
or in another way
Code:
   D(A_i)_new + D(A_b)_new = max{ D(A_i) + D(A_b)_new, sumOff }
which is an approximation of
Code:
    D(A)_new = max{D(A), sumOff},
update rhs with
Code:
    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
Code:
     D = D(A_i) + |D(A_b)|    ... (*)
relax diag
Code:
     D_new = max{D, sumOff} / alpha
update interal diag
Code:
     D(A_i)_new = D_new - D(A_b)  ...(**)
update rhs with
Code:
     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.
lengweee is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[Other] dynamicTopoFVMesh and pointDisplacement RandomUser OpenFOAM Meshing & Mesh Conversion 6 April 26, 2018 08:30
Problem with cyclic boundaries in Openfoam 1.5 fs82 OpenFOAM 36 January 7, 2015 01:31
[Commercial meshers] Fluent msh and cyclic boundary cfdengineering OpenFOAM Meshing & Mesh Conversion 48 January 25, 2013 04:28
Cyclic Boundary Condition Luiz Eduardo Bittencourt Sampaio (Sampaio) OpenFOAM Running, Solving & CFD 36 July 2, 2012 13:23
[mesh manipulation] Using createPatch in place of couplePatches sripplinger OpenFOAM Meshing & Mesh Conversion 8 November 13, 2009 08:14


All times are GMT -4. The time now is 09:48.