
[Sponsors] 
fvMatrix lduAddressing and coefficients calculation and accessing 

LinkBack  Thread Tools  Search this Thread  Display Modes 
September 8, 2016, 02:52 
fvMatrix lduAddressing and coefficients calculation and accessing

#1 
Member
Jerry
Join Date: Oct 2013
Location: Salt Lake City, UT, USA
Posts: 52
Rep Power: 12 
Dear Foamers,
This question has been raised for several times at least here in this Forum. It seems that they all didn't give a clear answer. That's why I am trying to bring it up again. http://www.cfdonline.com/Forums/ope...lculation.html http://www.cfdonline.com/Forums/ope...ddtscheme.html (The #7) http://www.cfdonline.com/Forums/ope...operator.html The calculation of the diagonal coefficients lots of times are done by the negSumDiag() function. The code is like this: Code:
void Foam::lduMatrix::negSumDiag() { const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower(); const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper(); scalarField& Diag = diag(); const labelUList& l = lduAddr().lowerAddr(); const labelUList& u = lduAddr().upperAddr(); for (register label face=0; face<l.size(); face++) { Diag[l[face]] = Lower[face]; Diag[u[face]] = Upper[face]; } } Code:
for (register label face=0; face<l.size(); face++) { Diag[l[face]] = Upper[face]; Diag[u[face]] = Lower[face]; } In santiagomarquezd's post, a 2D mesh was proposed: (Details were given in his post) *****  0  1  2  3  *****  4  5  6  7  ***** The calculation of diagonal coefficient for cell #1 is given in loops 0, 2 and 3 in the negSumDiag() function. Loop #0: Diag[u[0]] = Upper[0]; (Equivalent to: Diag[1] =Upper[0] ) Loop #2: Diag[l[2]] = Lower[2]; (Equivalent to: Diag[1] =Lower[2] ) Loop #3: Diag[l[3]] = Lower[3]; (Equivalent to: Diag[1] =Lower[3] ) In total, Diag[1] = (Upper[0] + Lower[2] + Lower[3]) This confuses me so much as it is very obvious that Diag[1] should be calculated as: Diag[1] = (Lower[0] + Upper[2] + Upper[3]) It might work for a symmetric matrix, but I really doubt about the nonsymmetric case. (correct me if I am wrong.) The other thing that makes me even more confusing is the lduMatrix::H function. The code: Code:
template<class Type> Foam::tmp<Foam::Field<Type> > Foam::lduMatrix::H(const Field<Type>& psi) const { tmp<Field<Type> > tHpsi ( new Field<Type>(lduAddr().size(), pTraits<Type>::zero) ); if (lowerPtr_  upperPtr_) { Field<Type> & Hpsi = tHpsi(); Type* __restrict__ HpsiPtr = Hpsi.begin(); const Type* __restrict__ psiPtr = psi.begin(); const label* __restrict__ uPtr = lduAddr().upperAddr().begin(); const label* __restrict__ lPtr = lduAddr().lowerAddr().begin(); const scalar* __restrict__ lowerPtr = lower().begin(); const scalar* __restrict__ upperPtr = upper().begin(); register const label nFaces = upper().size(); for (register label face=0; face<nFaces; face++) { HpsiPtr[uPtr[face]] = lowerPtr[face]*psiPtr[lPtr[face]]; HpsiPtr[lPtr[face]] = upperPtr[face]*psiPtr[uPtr[face]]; } } return tHpsi; } The part marked bold above means: (forget about the negtive sign in the code above and the source term below) we see that it's not using: Code:
for (register label face=0; face<nFaces; face++) { HpsiPtr[lPtr[face]] = lowerPtr[face]*psiPtr[lPtr[face]]; HpsiPtr[uPtr[face]] = upperPtr[face]*psiPtr[uPtr[face]]; } } Or there is some tricky part I didn't catch that makes these happen naturally. Discussions on this are really needed. Last edited by Jerryfan; September 8, 2016 at 15:32. 

September 8, 2016, 04:27 

#2  
Super Moderator
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51 
Dear Jerry,
first, good starting thread with all information and you took your time to prepare you question. I highly appreciate it. As you wrote me a private message, I will try to give you an answer based on the code that I see because I was not dealing with that classes till now. But first of all. If you think that the calculation is wrong, then all the validation cases that were calculated with openfoam should be wrong  luckily it is not like that. Okay lets start. You are asking about the behavior of the code of the negSumDiag() function. Can you please check you post because the code you posted and the code that is already implemented are identical. The second question belongs to the 2D case. Please check what you wrote again because the correction you mentioned is again similar to the one you mentioned above. Maybe you meant that one: Quote:
Code:
Diag[1] = (Lower[0] + Upper[2] + Upper[5]) Code:
face of Cell 1 has the faces (5 8 9 10) // Arbitrary numbering Code:
for (label face=0; face<nFaces; face++) { HpsiPtr[uPtr[face]] = lowerPtr[face]*psiPtr[lPtr[face]]; HpsiPtr[lPtr[face]] = upperPtr[face]*psiPtr[uPtr[face]]; } Code:
face of Cell 1 has the faces (5 8 9 10) // Arbitrary numbering fluxes of Cell 1 based on face ordering (1.243 12.353 1.529 8.982) I think the drawing you made is somehow a bit misleading because here we would use the ID of the faces either than moving them forward with pointers and addresses. In addition I tried to output the value of nFaces. Unfortunately the solvers icoFoam, simpleFoam and pimpleFoam are not calling this function so I do not know it exactly but as far as I understand, the variable face is an iterator through the flux list based on pointer and memory. Maybe it is not 100% correct what I told but it should be in that way... or if I would investigate more time, I will see that it is complete different. Due to the fact that there is no time, I hope that it is correct and help you. Furthermore based on the example you gave, the lower values should correspond to 0, 5 and the upper to 3 with respect to the cell 1. The last thing you mentioned may be correct but I do not know the correct equation. So I can not tell you anything.
__________________
Keep foaming, Tobias Holzmann 

September 8, 2016, 23:28 

#3  
Member
Jerry
Join Date: Oct 2013
Location: Salt Lake City, UT, USA
Posts: 52
Rep Power: 12 
Dear Tobi,
Thanks for your reply. I did make a mistake while I was copying and pasting. Sorry about this. This has been revised now. What I am trying to say is like this: Quote:
The loop inside the negSumDiag() function is controlled by face id rather than cell id, and the discretization coefficients contained in those Lower and Upper arrays are indexed by face id as well. As I mentioned above, Diag[1] will be visited three times. They are when: face=0, then u[0]=1 so that we have Diag[1], at the same time, Diag[1]=Upper[0]; face=2, then l[2]=1 so that we have a second Diag[1], at the same time, Diag[1]=Lower[2]; face=3, then l[3]=1 so that we have the third Diag[1], at the same time, Diag[1]=Lower[3]. That's why in total, Diag[1]=(Upper[0]+Lower[2]+Lower[3]) . But from the figure I attached above, it can be seen that Diag[1]=(Lower[0]+Upper[2]+Upper[3]) . why? I attached another figure, which is the matrix of the discretization coefficients. For cell #1, it's discretized equation will be: A10* T0 + A11* T1 + A12* T2 + A15 * T5 = b1 (For convenience, I just use T to represent the scalar variable we are solving). In the matrix, A10 corresponds to Lower[0], A12 corresponds to Upper[2] and A15 corresponds to Upper[3]. This can be found in the OpenFOAM guide/matrix. https://openfoamwiki.net/index.php/O...es_in_OpenFOAM However, if we use the code below in the negSumDiag() function, Quote:
Last edited by Jerryfan; September 9, 2016 at 01:09. Reason: the matrix image has been changed to a better one 

September 8, 2016, 23:37 

#4 
Member
Jerry
Join Date: Oct 2013
Location: Salt Lake City, UT, USA
Posts: 52
Rep Power: 12 
Too bad, the figures are not showing here. I don't know if you can see by right click and open it in another tab. If you cann't see it this way, I will think about other ways.


September 9, 2016, 00:04 

#5 
Member
Jerry
Join Date: Oct 2013
Location: Salt Lake City, UT, USA
Posts: 52
Rep Power: 12 
The images are showing now.


September 9, 2016, 00:32 

#6  
Member
Jerry
Join Date: Oct 2013
Location: Salt Lake City, UT, USA
Posts: 52
Rep Power: 12 
Dear Tobi,
It takes me quite a long time to prepare the above figures. Now let me explain what I mean for the lduMatrix::H function. First of all, the array lowerPtr() points to is not the flux. It is still the discretization coefficients as I explained above. I have to say the Lower() and Upper() coefficients in OpenFOAM (offdiagonal coefficients) are calculated correctly. I've followed the source code, they are done in a very appropriate way. That's why you get the correct "flux" (it's in fact, offdiagonal coefficients). I only have some doubts about the diagonal coefficients calculation as I explained above. What I am trying to point out here in the lduMatrix::H function is the indexes used for the HpsiPtr array. If we write out the corresponding mathematic equations of the negSumDiag() and the code snippet in the lduMatrix::H (forget about the source terms for H now), they are equivalent to: Let's still take the cell #1 in the 2D mesh above as an example. Quote:
face=0, then uPtr[0]=1 so that we have HpsiPtr[1], at the same time, HpsiPtr[1]=lowerPtr[0]*T[0]=A10*T[0]; face=2, then lPtr[2]=1 so that we have a second HpsiPtr[1], at the same time, HpsiPtr[1]=upperPtr[2]*T[2]=A12*T[2]; face=3, then lPtr[3]=1 so that we have the third HpsiPtr[1], at the same time, HpsiPtr[1]=upperPtr[3]*T[3]=A15*T[5]. In total, H[1] = (A10*T[0]+ A12*T[2]+ A15*T[5]). which is exactly what the equation represents. However, if we follow the way how the negSumDiag() function handles the indexes for the Diag array, we should write the code this way: Quote:
Last edited by Jerryfan; September 9, 2016 at 01:11. Reason: Equations are added 

September 11, 2016, 05:53 

#7 
Super Moderator
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,708
Blog Entries: 6
Rep Power: 51 
Hi,
I think I am not the right guy for your question. If I have time, I will investigate into that but up to now no idea. Do you have some literature in addition? maybe I will check Ferziger. By the way. When is th H function be called? as I checked out on Thursday, I did not figure out when this function is called. Kind regards, tobi
__________________
Keep foaming, Tobias Holzmann 

September 11, 2016, 08:44 

#8  
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 
Quote:
Quote:
The cite from Moukalled, Mangani, Darwish  The Finite Volume Method in Computational Fluid Dynamics An Advanced Introduction with OpenFOAMŪ and Matlab (Fluid Mechanics and Its Applications)  2015; pages 201202: Quote:
Quote:
It might be a thing related to the Hoperator used in SIMPLE algorithm 

September 23, 2016, 09:57 

#9  
New Member
Haixuan Ye
Join Date: May 2012
Location: Ann Arbor, MI, USA
Posts: 4
Rep Power: 14 
Quote:
But I guess I may misunderstand something, since all the places I can find these eqns write them in this way. 

January 4, 2017, 09:45 

#10  
New Member
Mohamed Elshahat Ouda
Join Date: May 2010
Posts: 29
Rep Power: 16 
Quote:
In case of advectiondiffusion equation the negSumDiag function use this smart trick to avoid calculating the diagonal coefficient again, alternatively it assemble it from the offdiagonal coefficients of the neighboring cells. Hope that helps. Kind regards. 

January 4, 2017, 16:21 

#11 
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 
Why? void Foam::lduMatrix::negSumDiag() is supposed to do a negative summation of the offdiagonal coefficients for each diagonal coefficient (each row in a matrix), isn't it? For each row you should only operate with coefficients in the same row.


January 4, 2017, 20:35 

#12  
New Member
Mohamed Elshahat Ouda
Join Date: May 2010
Posts: 29
Rep Power: 16 
Quote:
Code:
scalarField& Diag = diag(); Again this function do this smart trick to get the correct offdiagonal coefficient, namely coefficient from the corresponding neighbouring cell. Please correct me if this isn't true. Best regards. 

January 5, 2017, 16:44 

#13 
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 
By definition, for a specific diagonal coefficient A_ii, coefficients from the corresponding neighbouring cells are A_ij, not A_ji.


September 14, 2018, 18:55 

#14 
New Member
Chen Shen
Join Date: Sep 2018
Posts: 14
Rep Power: 7 
So I have written several lines of codes to construct a matrix and then call negSumDiag() to illustrate this problem and our confusion.
My mesh has 8 elements: int de; // for pause the program so I can read the output. scalarField& lower = UEqn.lower (); scalarField& diag = UEqn.diag (); scalarField& upper = UEqn.upper (); forAll (diag,i){ diag[i]=0; } forAll (upper,i){ upper[i]=i+8; } forAll (lower,i){ lower[i]=i+20; } Info << "diag \n" << diag << "upper\n" << upper << "lower" << lower; const lduAddressing & ldu = UEqn.lduAddr (); const labelUList& lo = ldu.lowerAddr (); const labelUList& up = ldu.upperAddr (); Info<<"lo:\n"<<lo<<"up:\n"<<up; scanf ("%d", &de); // input a int to go on. UEqn.negSumDiag(); Info << "diag \n" << diag << "upper\n" << upper << "lower" << lower; scanf ("%d", &de); ************************************************** ********************output************************ * diag 8{0} // diag is all 0. upper 12 ( 8 9 10 11 12 13 14 15 16 17 18 19 ) lower 12 ( 20 21 22 23 24 25 26 27 28 29 30 31 ) lo: 12 ( 0 0 0 1 1 2 2 3 4 4 5 6 ) up: 12 ( 1 2 4 3 5 3 6 7 5 6 7 7 ) So after setting values to the matrix stored in UEqn. Before UEqn.negSumDiag(); it is: 0.0 8.0 9.0 0 10.0 0 0 0 20.0 0.0 0 11.0 0 12.0 0 0 21.0 0 0.0 13.0 0 0 14.0 0 0 23.0 25.0 0.0 0 0 0 15.0 22.0 0 0 0 0.0 16.0 17.0 0 0 24.0 0 0 28.0 0.0 0 18.0 0 0 26.0 0 29.0 0 0.0 19.0 0 0 0 27.0 0 30.0 31.0 0.0 Lets look at the matrix values after UEqn.negSumDiag(); It is : 63.0 8.0 9.0 0 10.0 0 0 0 20.0 55.0 0 11.0 0 12.0 0 0 21.0 0 60.0 13.0 0 0 14.0 0 0 23.0 25.0 51.0 0 0 0 15.0 22.0 0 0 0 67.0 16.0 17.0 0 0 24.0 0 0 28.0 58.0 0 18.0 0 0 26.0 0 29.0 0 62.0 19.0 0 0 0 27.0 0 30.0 31.0 52.0 63=202122; 55=82324; and so on; I know the matrix is correctly printed because I am using gdbOF to print it. ( pfvMatrixFull UEqn f1, handy tool it you have not tried yet.) So my confusion is : UEqn.negSumDiag() sums the value in the same column for all diag components but not rows. This confuses me a lot. I think it should sum for the row because each row describe the relationship of each cell and its neighbors, not column, Tell me if the above is not clear and I will make it clear. Lets try to finish and close this question. Chen Last edited by Giantsda; September 17, 2018 at 18:02. 

September 28, 2018, 17:36 

#15 
New Member
Chen Shen
Join Date: Sep 2018
Posts: 14
Rep Power: 7 
I think I have figured this out:
1. void Foam::lduMatrix::negSumDiag() is negative sum its column , not its rows and then put the value to each diagonal entry. This can be confirmed by adding: fvScalarMatrix div (fvm::div (phi, T)); to almost any case and print the matrix of div and view it. The value of each diagonal entry is its negsum of its column. 2. Now for the question why we should sum for columns, not for rows, I did a simple case. The index is 1 based. and the matrix A is: Chen 

Tags 
fvmatrix, fvmatrix::h(), lduaddressing 
Thread Tools  Search this Thread 
Display Modes  

