[petsc-users] DMDA and ksp(ex46.c & ex22f.F90)

Smith, Barry F. bsmith at mcs.anl.gov
Mon Mar 11 23:24:54 CDT 2019



> On Mar 11, 2019, at 11:11 PM, Maahi Talukder <maahi.buet at gmail.com> wrote:
> 
> Hi
> Thank you so much for explanation.
> 
> So when you say connecting two points on the grid in the matrix entry, what do you mean by that? Do you mean that  matrix entry is calculated using two points ( (i,j) and (i',j') ) on the grid?

   I mean the value of the finite difference "stencil" that connects the two points.

   Take a look at src/ksp/ksp/examples/tutorials/ex22f.F90. The inner most loop

do 30,i=xs,xs+xm-1
          row(MatStencil_i) = i
          row(MatStencil_j) = j
          row(MatStencil_k) = k
     ....

sets one row of the matrix which has 7 nonzero columns corresponding to the 6 neighbors of the (i,j,k) point plus the diagonal entry which connects with itself.

>  And when you talked about thinking about the vector, how is it that the entry in the vector is a function of both i and j?

   For each (i,j) point on the mesh there is a vector entry. When indexing into vectors we use a single index value m; that has to be computed based on the two values of the logical coordinates on the mesh i and j.


>  Will it be the same if my dof = 1 ie solving for only one thing at each grid point?
   
   Yes it doesn't have anything to do with dof. The stencil_c value is used for multiple degrees of freedom per grid point.
> 
> 
> Regards,
> Maahi Talukder
> 
> On Mon, Mar 11, 2019 at 11:11 PM Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
> 
> 
> > On Mar 11, 2019, at 10:01 PM, Maahi Talukder <maahi.buet at gmail.com> wrote:
> > 
> > Hi
> > Thank you for your explanation. 
> > So is it so that it always connects two points? 
> 
>    Entries in matrices represent the connection between points in a vector (including the diagonal entries that are connections between a point and itself). 
> 
>    I think you are "over-thinking" MatStencil. It is just a handy way of managing the mapping from a grid point to an entry in a vector. And similarly from two grid points and an entry in a matrix.
> 
> > I am still hazy on the concept. Could you suggest me some more elaborate book or something to go through so that I can grasp the whole concept? How about your book ' Domain Decomposition'? Does it cover this topic?
> 
>   No, it doesn't discuss this issue.
> 
>    Barry
> 
> > Please let me know
> > 
> > Regards,
> > Maahi Talukder 
> > 
> > On Mon, Mar 11, 2019 at 9:46 PM Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
> > 
> > 
> > > On Mar 11, 2019, at 7:07 PM, Maahi Talukder via petsc-users <petsc-users at mcs.anl.gov> wrote:
> > > 
> > > 
> > > Thank you for your reply. 
> > > 
> > > I still have some confusion. So if (i,j) is a point on the structured grid( Where "i" is the column and "j" is the row), and the information associated with the (i,j) point on the grid is stored in some (m,n) location of the matrix A (Where Ax =b),
> > 
> >    Almost. Think about the vector (not the matrix) first. The point (i,j) on the mesh has location m in the vector (m is a function of both i and j) Now think about another point on the mesh (say next to the first point), call it (i',j') this point also has a location in the vector, call it n. Now consider the matrix entry containing the value that connects m and n. On the mesh it is associated with the point (i,j) AND the point (i',j') hence the row in the matrix is a function of the first point (i,j) while the column is a function of the second point (i',j'). 
> > 
> >   Now the diagonal entry in matrix a_{mm} has (i,j) for both "row" and "column" stencils, but the off-diagonal entries a_{mn} has (i,j) for the row but (i',j') for the column. 
> > 
> > > I still don't 
> > > understand why both of  row(MatStencil_i,1) and row(MatStencil_j,1) are necessary? I mean is it something like mapping "i" from grid to its location in the matrix?
> > 
> > 
> > > Would you please explain that?
> > 
> > 
> > > 
> > > Regards,
> > > Maahi 
> > > 
> > > On Mon, Mar 11, 2019 at 4:41 PM Patrick Sanan <patrick.sanan at gmail.com> wrote:
> > > There are two different types of rows and columns:
> > > 1. Rows and columns in a grid
> > > 2. Rows and columns in a matrix
> > > 
> > > "i" and "j"  refer to rows and columns in the grid, but "row" and "col"  refer to rows and columns in the matrix.
> > > 
> > > 
> > > 
> > > Am Mo., 11. März 2019 um 21:18 Uhr schrieb Maahi Talukder via petsc-users <petsc-users at mcs.anl.gov>:
> > > Hello all, 
> > > 
> > > I am trying to solve Poisson Equation on structured grid using 9-point stencil in 2D. Now to setup my matrix, I came across C structure MatStencil in ex22f.F90
> > > 
> > > ...........................................................................................................
> > > call DMDAGetCorners
> > > (da,xs,ys,zs,xm,ym,zm,ierr)
> > > 
> > > 
> > > 107:       do
> > >  10,k=zs,zs+zm-1
> > > 
> > > 108:         do
> > >  20,j=ys,ys+ym-1
> > > 
> > > 109:           do
> > >  30,i=xs,xs+xm-1
> > > 
> > > 110: 
> > >           row(MatStencil_i) = i
> > > 
> > > 111: 
> > >           row(MatStencil_j) = j
> > > 
> > > 112: 
> > >           row(MatStencil_k) = k
> > > 
> > > 113:           if
> > >  (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1 .or. k.eq.mz-1) then
> > > 
> > > 114: 
> > >             v(1) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
> > > 
> > > 115:             call MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES
> > > ,ierr)
> > > 
> > > 116:           else
> > > 117: 
> > >             v(1) = -HxHydHz
> > > 
> > > 118: 
> > >              col(MatStencil_i,1) = i
> > > 
> > > 119: 
> > >              col(MatStencil_j,1) = j
> > > 
> > > 120: 
> > >              col(MatStencil_k,1) = k-1
> > > 
> > > 121: 
> > >             v(2) = -HxHzdHy
> > > 
> > > 122: 
> > >              col(MatStencil_i,2) = i
> > > 
> > > 123: 
> > >              col(MatStencil_j,2) = j-1
> > > 
> > > 124: 
> > >              col(MatStencil_k,2) = k
> > > 
> > > 125: 
> > >             v(3) = -HyHzdHx
> > > 
> > > 126: 
> > >              col(MatStencil_i,3) = i-1
> > > 
> > > 127: 
> > >              col(MatStencil_j,3) = j
> > > 
> > > 128: 
> > >              col(MatStencil_k,3) = k
> > > 
> > > 129: 
> > >             v(4) = 2.0*(HxHydHz + HxHzdHy + HyHzdHx)
> > > 
> > > 130: 
> > >              col(MatStencil_i,4) = i
> > > 
> > > 131: 
> > >              col(MatStencil_j,4) = j
> > > 
> > > 132: 
> > >              col(MatStencil_k,4) = k
> > > 
> > > 133: 
> > >             v(5) = -HyHzdHx
> > > 
> > > 134: 
> > >              col(MatStencil_i,5) = i+1
> > > 
> > > 135: 
> > >              col(MatStencil_j,5) = j
> > > 
> > > 136: 
> > >              col(MatStencil_k,5) = k
> > > 
> > > 137: 
> > >             v(6) = -HxHzdHy
> > > 
> > > 138: 
> > >              col(MatStencil_i,6) = i
> > > 
> > > 139: 
> > >              col(MatStencil_j,6) = j+1
> > > 
> > > 140: 
> > >              col(MatStencil_k,6) = k
> > > 
> > > 141: 
> > >             v(7) = -HxHydHz
> > > 
> > > 142: 
> > >              col(MatStencil_i,7) = i
> > > 
> > > 143: 
> > >              col(MatStencil_j,7) = j
> > > 
> > > 144: 
> > >              col(MatStencil_k,7) = k+1
> > > 
> > > 145:       call MatSetValuesStencil(jac,i1,row,i7,col,v,INSERT_VALUES
> > > ,ierr)
> > > 
> > > 146:           endif
> > > .....................................................................................
> > > What I am confused about is what it means to have the value of row in i and j directions(row(MatStencil_i,1) & row(MatStencil_j,1)).
> > > Same confusion goes for the column values as well. I mean generally in a  2D Matrix row values are in j/y direction and column values are in i/x direction.
> > > Could you please explain that?
> > > 
> > > Regards,
> > > Maahi Talukder
> > > Department of Mechanical Engineering
> > > Clarkson University
> > 
> 



More information about the petsc-users mailing list