
[Sponsors] 
Information about equation relaxation (patch contribution) 

LinkBack  Thread Tools  Search this Thread  Display Modes 
April 6, 2017, 04:55 
Information about equation relaxation (patch contribution)

#1 
Super Moderator

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):
Hope some expert in that topic can tell me whats going on here. Code:
else { // For noncoupled 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]); } }
__________________
Keep foaming, Tobias Holzmann Last edited by Tobi; April 6, 2017 at 07:09. Reason: small mistake  we remove and add the boundary stuff to the diag 

April 6, 2017, 06:44 

#2 
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 863
Rep Power: 23 
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. 

April 6, 2017, 07:17 

#3 
Super Moderator

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 

April 6, 2017, 11:03 

#4 
Super Moderator

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 

April 6, 2017, 16:06 

#5 
Senior Member

Hi,
Different times are the result of these lines (https://cpp.openfoam.org/v4/a00885.h...f646011c212e): Code:
if (alpha <= 0) { return; } 

April 6, 2017, 16:27 

#6 
Super Moderator

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 } b) The matrix represents a nonscalar matrix and the boundary values in the direction differs. Why? Because first we add the magnitue of the maximum component: Code:
558 // For noncoupled 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 } Code:
648 forAll(pa, face) 649 { 650 D[pa[face]] = cmptMin(iCoeffs[face]); 651 } 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 

April 6, 2017, 19:05 

#7 
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 863
Rep Power: 23 
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. 

August 25, 2018, 04:28 

#8 
New Member
wei leng
Join Date: Mar 2014
Posts: 3
Rep Power: 8 
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.cfdonline.com/Forums/openfoam/167028extractingmatrixdatafvvectormatrix.html#post703887 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 Code:
D(A)_new = max{D(A), sumOff} / alpha, Code:
f = f + (D(A)_new  D(A)) * uold 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, Code:
f_b,new = f_b + ( D(A_b)_new  D(A_b) ) * x update diag with Code:
D(A_i)_new = max{ D(A_i), sumOff  D(A_b)_new } Code:
D(A_i)_new + D(A_b)_new = max{ D(A_i) + D(A_b)_new, sumOff } Code:
D(A)_new = max{D(A), sumOff}, 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) ... (*) Code:
D_new = max{D, sumOff} / alpha Code:
D(A_i)_new = D_new  D(A_b) ...(**) Code:
f_i,new = f_i + ( D(A_i)_new  D(A_i) ) * x 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. 

Thread Tools  Search this Thread 
Display Modes  


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