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

Matrix storage and diagonal coeffients calculation

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

Like Tree8Likes
  • 2 Post By hjasak
  • 1 Post By Tobermory
  • 1 Post By Tobermory
  • 1 Post By nipinl
  • 1 Post By nipinl
  • 1 Post By Rvadrabade
  • 1 Post By nipinl

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 5, 2010, 13:40
Default Matrix storage and diagonal coeffients calculation
  #1
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 452
Rep Power: 22
santiagomarquezd will become famous soon enough
Hi all, here we're dealing with lduMatrix coefficients storage and calculation. In classical books for FVM like the one written by Patankar the diagonal coefficient arises as the summation of all the neighbour's coefficients, i.e. the coefficients in the same row, talking about matrices. But looking inside in the implementation of laplacian and full upwind schemes, both uses the negSumDiag method, that sums by column not by row to obtain the diagonal coefficient. Might be the matrices transposed for storage?

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Research Scientist
Research Center for Computational Methods (CIMEC) - CONICET/UNL
Tel: 54-342-4511594 Int. 7032
Colectora Ruta Nac. 168 / Paraje El Pozo
(3000) Santa Fe - Argentina.
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   February 6, 2010, 07:31
Default
  #2
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,881
Rep Power: 31
hjasak will become famous soon enough
Nope: the negSumDiag() will sum up on rows. What may be confusing you is the names of the addressing arrays, but the operation is correct:
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 unallocLabelList& l = lduAddr().lowerAddr();
const unallocLabelList& u = lduAddr().upperAddr();

for (register label face=0; face<l.size(); face++)
{
Diag[l[face]] -= Lower[face];
Diag[u[face]] -= Upper[face];
}
}

Lower stand for "lower addressing index"; since the addressing is done through the upper triangle, the lower index is the index for a row.

Enjoy,

Hrv
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   February 8, 2010, 10:15
Default
  #3
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 452
Rep Power: 22
santiagomarquezd will become famous soon enough
I know I'm wrong, but I can't figure where is my error, I can't completely understand the phrase "Lower stand for "lower addressing index"; since the addressing is done through the upper triangle, the lower index is the index for a row." Lower is Upper Triangular? Maybe an example can help my point of view and my error. Suppose we have a mesh with two layers of four hexahedron, representing a 2D mesh.

*-------*-------*-------*-------*
| 0 | 1 | 2 | 3 |
*-------*-------*-------*-------*
| 4 | 5 | 6 | 7 |
*-------*-------*-------*-------*

the sides of this hexahedron have unitary areas. When the mesh is obtained we have the following arrays:

Owner={0, 0, 1, 1, 2, 2, 3, 4, 5, 6, ...}
Neigbour={1, 4, 2, 5, 3, 6, 7, 5, 6, 7}

Using an advective velocity of 1 in horizontal direction (to the right) and applying the FVM method we can assemble for cell 1 (equation #1, starting in 0) the following equation:

-phi_0 + phi_1=0

due the direction of the velocity, the negative coefficient (-1) will be ever in the lower triangular.
Running the example with scalarTransportFoam and debugging (in my case with gdb) we can can stop the program when negSumDiag is called within the advective term calculation.
At this point we have the following advective matrix (M):

0 0 0 0 0 0 0 0
-1 0 0 0 0 0 0 0
0 -1 0 0 0 0 0 0
0 0 -1 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 -1 0 0 0
0 0 0 0 0 -1 0 0
0 0 0 0 0 0 -1 0

l array has: {0, 0, 1, 1, 2, 2, 3, 4, 5, 6}, so l (from lowerAddr method) has the same data of Owner (at least the first 10 terms),
u array has: {1, 4, 2, 5, 3, 6, 7, 5, 6, 7}, so u (from upperAddr method) has the same data of Neighbour.
On the other hand, we have for Lower: {-1, -0, -1, -0, -1, -0, -0, -1, -1, -1} (from the lower attribute of the class), it says that Lower is the Lower Triangular, and for Upper: {0, 0, 0, 0, 0, 0, 0, 0, 0, 0} (from the upper attribute of the class), the Upper Triangular. Finally for Diag: {0, 0, 0, 0, 0, 0, 0, 0}.
Looking at the Lower Triangular l array give us the columns of elements, and u the rows. So for element number 0, i.e. -1, we have row: 1, column: 0, like in the matrix.
Running the loop in negSumDiag, for example for the equation presented before (equation #1), calculation of diagonal coefficient is given in loops 0, 2 and 3 (i.e. face index in loop takes values 0, 2 and 3), in loop 0 we have:

Diag[l[0]] -= Lower[0];
Diag[u[0]] -= Upper[0];

Second line operates on Diag[1] M[1][1] subtracting the first element of Upper, i.e M[0][1] (row is given by l and column by u). In loop 2:

Diag[l[2]] -= Lower[2];
Diag[u[2]] -= Upper[2];

First line operates on M[1][1] subtracting M[2][1] (the third element of Lower whose row is given by u and column by l). Finally in loop 3 we have M[1][1]=M[1][1]-M[5][1]. Finally the diagonal coefficient 1 was calculated from elements in the same column but in different rows when I had expected that this calculation would be made from elements in the same row. i.e. M[1][1]=M[1][1]-M[1][0].

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Research Scientist
Research Center for Computational Methods (CIMEC) - CONICET/UNL
Tel: 54-342-4511594 Int. 7032
Colectora Ruta Nac. 168 / Paraje El Pozo
(3000) Santa Fe - Argentina.
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   February 11, 2010, 13:00
Default
  #4
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 452
Rep Power: 22
santiagomarquezd will become famous soon enough
Any clues...? Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Research Scientist
Research Center for Computational Methods (CIMEC) - CONICET/UNL
Tel: 54-342-4511594 Int. 7032
Colectora Ruta Nac. 168 / Paraje El Pozo
(3000) Santa Fe - Argentina.
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   March 26, 2010, 07:28
Default
  #5
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 452
Rep Power: 22
santiagomarquezd will become famous soon enough
We've studying the matter for an extra while and found that row summation is not done in the same way that in Patankar's book, but the results are correct. Column summation that I explained in past posts is only an algorithmic implementation for the sake of simplification in machine use. Explanation is large, if somebody needs more details, contact me or contact nisi as well. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Research Scientist
Research Center for Computational Methods (CIMEC) - CONICET/UNL
Tel: 54-342-4511594 Int. 7032
Colectora Ruta Nac. 168 / Paraje El Pozo
(3000) Santa Fe - Argentina.
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   November 3, 2014, 14:07
Default Could you please post the explanation
  #6
Member
 
Ganesh Vijayakumar
Join Date: Jan 2010
Posts: 44
Rep Power: 14
ganeshv is on a distinguished road
Hi!

I'm so thankful for finding this post. I've been wondering about this for a while. I see that this implementation is ok for symmetric matrices. But I still don't see it for asymmetric matrices. Could you please post the explanation for all of us.

I also have a follow up question. Even if this implementation works for symmetric matrices... why would it help so much to implement it this way? What do we save by introducing this (for lack of a better word) "confusion" ?

ganesh
ganeshv is offline   Reply With Quote

Old   March 23, 2020, 19:17
Default
  #7
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 797
Rep Power: 15
sharonyue is on a distinguished road
This topic was discussed at length in this link: http://www.cfd-china.com/topic/547
in Chinese though...
__________________
My OpenFOAM algorithm website: http://dyfluid.com
By far the largest Chinese CFD-based forum: http://www.cfd-china.com/category/6/openfoam
We provide lots of clusters to Chinese customers, and we are considering to do business overseas: http://dyfluid.com/DMCmodel.html
sharonyue is offline   Reply With Quote

Old   May 23, 2020, 06:09
Default
  #8
Member
 
Rahul Vadrabade
Join Date: Apr 2018
Posts: 41
Rep Power: 6
Rvadrabade is on a distinguished road
Few comments.

As everyone is aware of equation after discretization A_{P} \psi_{P}+\sum_{N} A_{N} \psi_{N}=\mathbf{S}_{p}.

For div(phi,U)

A_{P}=-\sum_{N} A_{N}+totalFlux

for 2d: A_{P}=A_{W}+A_{E}+A_{S}+A_{N}+\{phi_e-phi_w +phi_n-phi_s\}

For linear (w=0.5) and upwind (w=1) {w=weight}

A_{W} = -w\bullet phi_{w}
A_{E} = (1-w)\bullet phi_{e}
A_{N} = (1-w)\bullet phi_{n}
A_{S} = -w\bullet phi_{s}

So, yes it seems that OF does sum by column but it is not.


Reference:- Section 5.7.2 Hybrid differencing scheme for multi-dimensional convection-diffusion;
An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Book by H K Versteeg and W Malalasekera

Note that:- Coefficient in the book are calculated assuming equation A_{p} \psi_{p}=\sum_{N} A_{N} \psi_{N}+\mathbf{S}_{p} so they have opposite sign but taken care at matrix construction.
Rvadrabade is offline   Reply With Quote

Old   May 23, 2020, 12:22
Default
  #9
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,881
Rep Power: 31
hjasak will become famous soon enough
First, you’re wrong because FOAM does this in face addressing so you don’t have w and (1-w)

Second, this has worked for 27 years without a change, producing perfect results. Did you actually produce the test case that reproduces the error?
Ardali and Rvadrabade like this.
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   May 24, 2020, 01:30
Default
  #10
Member
 
Rahul Vadrabade
Join Date: Apr 2018
Posts: 41
Rep Power: 6
Rvadrabade is on a distinguished road
Hi Prof. Hrvoje Jasak

I agree with you and OF implementation is absolutely correct, i don't have any objections.

I would like to convey is that, the diagonal coefficient is not just addition of off diagonal coeffs but totalFlux term should be considered.

for 2D:
A_{P}=A_{W}+A_{E}+A_{S}+A_{N}+\{phi_e-phi_w +phi_n-phi_s\}

With the face addressing it is taken care however with cell addressing (as per books ) the extra term totalFlux needs attention.

Also, regarding use of 'w and (1-w)' probably one can refer here
https://www.openfoam.com/documentati...n-details.html

Thank you
Rvadrabade is offline   Reply With Quote

Old   December 23, 2020, 06:25
Default
  #11
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 456
Rep Power: 11
Tobermory will become famous soon enough
Quote:
Originally Posted by santiagomarquezd View Post
We've studying the matter for an extra while and found that row summation is not done in the same way that in Patankar's book, but the results are correct. Column summation that I explained in past posts is only an algorithmic implementation for the sake of simplification in machine use. Explanation is large, if somebody needs more details, contact me or contact nisi as well. Regards.
Dear Santiago - I have been trying to puzzle this out for a few days ... could you share your insight on how this works? Thanks & regards.
Rvadrabade likes this.
Tobermory is offline   Reply With Quote

Old   January 5, 2021, 04:49
Default
  #12
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 456
Rep Power: 11
Tobermory will become famous soon enough
Okay, I have finally puzzled this out - the coding in negSumDiag is in effect a piece of "sleight of hand" that makes for efficient coding.

Here is the key: when you write out the matrix coefficients for the divergence term, employing interpolation of the variable to the face and the face flux with the relevant sign, you get a set of off-diagonal and diagonal coefficients for the master and slave cell (owner and neighbour in OF terminology, although the term neighbour is a bit ambiguous) of each face. The trick is to realise that the diagonal coeffs are equal to the negative of the off-diagonal coeffs of the neighbouring cell, and not their own cell. This is why negSumDiag appears to "sum the columns in the matrix" ... it actually is!

I was confused for a while since I had written down the divergence in non-conservative form for some crazy reason, in which case the diagonal coeffs are then the negatives of their own cell's off-diagonal coefficients, and negSumDiag would then have had to sum the rows.
Rvadrabade likes this.

Last edited by Tobermory; January 5, 2021 at 11:49.
Tobermory is offline   Reply With Quote

Old   January 11, 2022, 13:15
Default confusion with diag term
  #13
New Member
 
Nipin L
Join Date: Nov 2012
Location: India
Posts: 14
Rep Power: 11
nipinl is on a distinguished road
I made a 2D laminar flow with 3 X 3 cells. Printing out the div(phi,U) as an fvMatrix produced the following output:
divUEqn:true true true 12{-0} 9(1e-05 0 0 1e-05 0 0 1e-05 0 0) 12{0}
Clearly, upper and lower are zeros, but diag is having non-zero values.


I used upwind scheme. The 1e-5 value is matching with flux (U at inlet = 0.001, Sf =0.01, phi = Sf*U = 1e-5. But my doubt is, where exactly the diag is getting values other than negSumDiag()?
Rvadrabade likes this.

Last edited by nipinl; January 14, 2022 at 04:51. Reason: corrected the expression for flux
nipinl is offline   Reply With Quote

Old   January 11, 2022, 22:06
Default Few references
  #14
Member
 
Rahul Vadrabade
Join Date: Apr 2018
Posts: 41
Rep Power: 6
Rvadrabade is on a distinguished road
It is formulated cleverly. I also find difficulties understanding the concept but following blogs gave me better understanding of the topics. I hope this helps if you are not aware of.

LaplacianFoam and fvmatrix -> https://zhuanlan.zhihu.com/p/32940120

IcoFoam data structure -> https://zhuanlan.zhihu.com/p/32679170

Few tools that might help you while debugging:-

https://github.com/Rvadrabade/Debugg...al-Studio-Code

gdbOF ->https://gitlab.com/flowcrunchpublic/...e/RahulIssues1
Rvadrabade is offline   Reply With Quote

Old   January 11, 2022, 23:50
Default Thank you
  #15
New Member
 
Nipin L
Join Date: Nov 2012
Location: India
Posts: 14
Rep Power: 11
nipinl is on a distinguished road
Thank you Rahul for your quick reply and the links for further read. Actually I was going through the gdbOF, is it compatible with v2012?
Rvadrabade likes this.
nipinl is offline   Reply With Quote

Old   January 12, 2022, 09:26
Default
  #16
Member
 
Rahul Vadrabade
Join Date: Apr 2018
Posts: 41
Rep Power: 6
Rvadrabade is on a distinguished road
Yes, ofcourse it is. Let me know if you face any issue.
nipinl likes this.
Rvadrabade is offline   Reply With Quote

Old   January 14, 2022, 04:48
Default confusion with diag term clarified
  #17
New Member
 
Nipin L
Join Date: Nov 2012
Location: India
Posts: 14
Rep Power: 11
nipinl is on a distinguished road
Quote:
Originally Posted by nipinl View Post
I made a 2D laminar flow with 3 X 3 cells. Printing out the div(phi,U) as an fvMatrix produced the following output:
divUEqn:true true true 12{-0} 9(1e-05 0 0 1e-05 0 0 1e-05 0 0) 12{0}
Clearly, upper and lower are zeros, but diag is having non-zero values.


I used upwind scheme. The 1e-5 value is matching with flux (U at inlet = 0.001, Sf =0.1phi = Sf*U = 1e-4. div(phi, U) = phi*U*Sf= 1e-5.). But my doubt is, where exactly the diag is getting values other than negSumDiag()?

Infact, I used bounded upwind. In bounded div schemes, an additional term is added to accelerate the convergence(explained here). boundedConvectionScheme.C at line 76 reads

template<class Type>

tmp<fvMatrix<Type>>

boundedConvectionScheme<Type>::fvmDiv

(

const surfaceScalarField& faceFlux,

const GeometricField<Type, fvPatchField, volMesh>& vf

) const

{

return

scheme_().fvmDiv(faceFlux, vf)

- fvm::Sp(fvc::surfaceIntegrate(faceFlux), vf);

}
In my case, the last line (surfaceIntegrate(faceFlux), vf)) added inlet fluxes to the diagonal. (looked like it replaced the diagonal entries for cells adjacent to inlet. I need to dig out further here)

When I changed the div scheme to Gauss linear, my diag() matched with negSumDiag() .



@Rahul, the links are super useful. Thank you for that :-)
Rvadrabade likes this.
nipinl 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
Diagonal storage for c-grid airfoil possible? zonexo Main CFD Forum 1 May 10, 2006 13:43


All times are GMT -4. The time now is 23:39.