# Matrix storage and diagonal coeffients calculation

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

 February 5, 2010, 14:40 Matrix storage and diagonal coeffients calculation #1 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 418 Rep Power: 15 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. Post-doctoral Fellow Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL T.E.: 54-342-4511594 Ext. 1005 Güemes 3450 - (3000) Santa Fe Santa Fe - Argentina http://www.cimec.org.ar

 February 6, 2010, 08:31 #2 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,763 Rep Power: 21 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(*this).lower(); const scalarField& Upper = const_cast(*this).upper(); scalarField& Diag = diag(); const unallocLabelList& l = lduAddr().lowerAddr(); const unallocLabelList& u = lduAddr().upperAddr(); for (register label face=0; face

 February 8, 2010, 11:15 #3 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 418 Rep Power: 15 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. Post-doctoral Fellow Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL T.E.: 54-342-4511594 Ext. 1005 Güemes 3450 - (3000) Santa Fe Santa Fe - Argentina http://www.cimec.org.ar

 February 11, 2010, 14:00 #4 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 418 Rep Power: 15 Any clues...? Regards. __________________ Santiago MÁRQUEZ DAMIÁN, Ph.D. Post-doctoral Fellow Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL T.E.: 54-342-4511594 Ext. 1005 Güemes 3450 - (3000) Santa Fe Santa Fe - Argentina http://www.cimec.org.ar

 March 26, 2010, 08:28 #5 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 418 Rep Power: 15 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. Post-doctoral Fellow Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL T.E.: 54-342-4511594 Ext. 1005 Güemes 3450 - (3000) Santa Fe Santa Fe - Argentina http://www.cimec.org.ar

 November 3, 2014, 15:07 Could you please post the explanation #6 Member   Ganesh Vijayakumar Join Date: Jan 2010 Posts: 43 Rep Power: 7 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

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post zonexo Main CFD Forum 1 May 10, 2006 13:43

All times are GMT -4. The time now is 03:22.