> On Sep 11, 2022, at 9:56 AM, Tu, Jiannan <jiannan...@uml.edu> wrote: > > Barry, > > Thank you. > > Yes, the DMDA is created with periodic BC. But it might not be the periodic > BC since in azimuthal direction, the ghost grid -1 is actually grid Nmax, and > Nmax+1 is grid 0.
Consider the following, one has n points from 0 to n-1 where n = 4, so I think with your definition of Nmax, Nmax is 3, Nmax + 1 is 4, -1 0 1 2 3 4 + * * * * + Thus there is one grid point at each end. The * represent regular grid points and the + ghost locations. Is this your situation? This is what we call periodic with a stencil width of 1 with DMDA and the DMCreateMatrix() should automatically provide the correct preallocation and nonzero structure. Can you please do a MatView() on the result from DMCreateMatrix() for a very small grid before you set any values and you should see the locations for the coupling across the periodic boundary. It sounds like you have a more complicated cross periodic coupling, we need to understand it exactly so we can determine how to proceed. Barry > > Still the problem is how to properly pre-allocate matrix. Otherwise either > “Augument out of range” error or code stuck at inserting values if > MatSetOption() is called. > > Jiannan > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Sunday, September 11, 2022 12:10 AM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Cc: petsc-users <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> > Subject: Re: [petsc-users] Using matrix-free with KSP > > CAUTION: This email was sent from outside the UMass Lowell network. > > > Have you created the DMDA with the periodic boundary conditions arguments? > Or are your boundary conditions not simply periodic? > > Compare DMCreateMatrix_DA_3d_MPIAIJ() and > DMCreateMatrix_DA_3d_MPIAIJ_Fill() in src/dm/impls/da/fdda.c I think they > both attempt to handle the periodic case. > > Barry > > > > > > On Sep 10, 2022, at 10:14 PM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > Barry, > > I think I understand how to set values for the intra-grid coupling. But I am > not clear about the way to set values for inter-grid coupling. When a > component couples to itself at above, below, right, left, front and back > grids, diagonal is set as 1 for the second array of DMDASetBlockFills(). How > about if involving more grids? How does DMDA know where to allocate memory? > In my case, a component couples to itself at the other end of the grid > because of periodic boundary condition. There is an error message “Augment > out of range . Insert new nonzero at global row/column (26, 520526) into > matrix” because of this far away coupling. The column index 520526 is related > to the other end of the grid in azimuthal direction. > > Thanks, > Jiannan > > > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Saturday, September 10, 2022 1:10 PM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Cc: petsc-users <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> > Subject: Re: [petsc-users] Using matrix-free with KSP > > CAUTION: This email was sent from outside the UMass Lowell network. > > > If you are using DMCreateMatrix() then MatMPIAIJSetPreallocation is not > needed, it is automatically handled by the DM. > > It sounds like that the values you set in DMDASetBlockFills() did not > include all the coupling, hence the preallocation was wrong, hence the time > to insert values was too long. > > Barry > > > > > On Sep 10, 2022, at 8:16 AM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > Just found MatMPIAIJSetPreallocation is needed. Now Jacobin insertion and > assembly is fast. > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Thursday, September 8, 2022 9:53 PM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Subject: Re: [petsc-users] Using matrix-free with KSP > > CAUTION: This email was sent from outside the UMass Lowell network. > > > > > > > On Sep 8, 2022, at 8:34 PM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > Barry, > > Thank you. > > I set up these two arrays according to the non-zero pattern of the Jacobian. > The DMDASetBlockFills() is then called after the call to DMCreateMatrix(da, > &A). > > Absolutely not! It must be called before DMCreateMatrix() so that the DM > knows how to construct exactly the matrix you need. > > > > > > Nevertheless, the program execution is still killed with exit code 9, which I > assume is due to overuse of the memory. My desktop machine has 32 GB RAM. > > I ran the code with a smaller mesh, 100x45x90, the result is the same. With > that number of grids, 26 unknowns, and up to 10 non-zero at each matrix row, > the estimated memory to store 100x45x90x26x10 elements should be less than 1 > GB (double precision). I am wondering why the matrix still takes too much > memory. Maybe I have to use matrix-free? > > Thank you, > Jiannan > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Thursday, September 8, 2022 4:19 PM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Subject: Re: [petsc-users] Using matrix-free with KSP > > CAUTION: This email was sent from outside the UMass Lowell network. > > > > > > > > On Sep 8, 2022, at 3:39 PM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > Barry, > > Thank you very much for the detailed description. > > So the first array, dfill, is for intra grids and the second one, ofill is > for the inter-grid coupling? > > In the case for inter-grid coupling, a component is only coupled to itself, > say x in your example is coupled to itself on above, below, right, left, > front, and back, how to set values for the second array? > > Then the off-diagonal one would just have values on its diagonal. > > > > > > Jiannan > > From: Barry Smith <mailto:bsm...@petsc.dev> > Sent: Thursday, September 8, 2022 2:12 PM > To: Tu, Jiannan <mailto:jiannan...@uml.edu> > Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> > Subject: Re: [petsc-users] Using matrix-free with KSP > > CAUTION: This email was sent from outside the UMass Lowell network. > > > > > > > > > On Sep 8, 2022, at 1:47 PM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > Barry, > > Thank you very much. > > DMDASetBlockFills() needs two input arrays specifying sparsity pattern. In my > case, I have 26 components at each grid. Is it correct that I need to have > two arrays of 26x26 elements with value 1 for coupling and 0 for no coupling? > > Exactly > > > > > > > Also how to determine when the value should be 1 or 0? > > That is problem dependent. Given paper-and-pencil you know for your stencil > which components are coupled with other components at each grid point and > which components are coupled with components at neighboring grid points. > > Note the second array contains information about all the neighboring grid > points. Above, below, right, left, front back. Normally the coupling is > symmetric so the nonzero structures of those 6 blocks are identical. If, for > you particular problem, the nonzero structures of the six are not identical > you use the union of the nonzero structures of the six of them. > > > 2d example say you have x,y,z at each grid point and the inter-point > coupling is x is connected to y and y is connect to x then the first array is > > [ 1 1 0 ; 1 1; 0 ; 0 0 1]; the zeros are because x is not coupled to z and y > is not coupled to z. > > For intra grid points, say x is coupled to z and itself, z is coupled to x > and itself and y is not coupled to itself then > > [1 0 1; 0 0 0; 1 0 1]; the middle row is zero because y is not coupled to > anything and the zeros in the first and last are because x and z are not > coupled to y. > > > > > > > > > Each component involves 7 stencils (i-1, j, k), (I, j, k), (i+1, j, k), (I, > j-1, k), (I, j+1, k), (I, j, k-1), and (I, j, k+1), and couples with several > other components. > > Jiannan > > > > Sent from Mail > <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgo.microsoft.com%2Ffwlink%2F%3FLinkId%3D550986&data=05%7C01%7CJiannan_Tu%40uml.edu%7C68957aad881840dcf87f08da93ab77f6%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637984661923188386%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=vWNBlC%2BrqWkX5Q0%2Bh6XVncZHaaqFAgShn2RthwcJem4%3D&reserved=0> > for Windows > > From: Barry Smith <mailto:bsm...@petsc.dev> > Sent: Wednesday, September 7, 2022 11:53 AM > To: Tu, Jiannan <mailto:jiannan...@uml.edu> > Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> > Subject: Re: [petsc-users] Using matrix-free with KSP > > CAUTION: This email was sent from outside the UMass Lowell network. > > > DMDASetBlockFills() or DMDASetBlockFillsSparse() (they are just different > ways of providing the same information) will help you here enormously, your > type of case is exactly what they were designed for. > > Barry > > > > > > > > > On Sep 7, 2022, at 10:47 AM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > Barry and Hong, > > Thank you. > > There are 26 components at each grid and there are not fully coupled in terms > of stiff functions. Mutual coupling is among about 6 components. > > I would prefer not using matrix-free since the Jacobina is not difficult to > calculate and only up to 10 non-zeros at each row. I'll try > DMDASetBlockFills() or DMDASetBlockFillsSparse() and see how they can reduce > the memory usage. > > Jiannan > <E6340AADF8494DB7861EA4659E7713CE.png> > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Tuesday, September 6, 2022 11:33 PM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> > <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> > Subject: Re: [petsc-users] Using matrix-free with KSP > > CAUTION: This email was sent from outside the UMass Lowell network. > > > > > > > > > > On Sep 6, 2022, at 11:00 PM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > I am using TS IMEX to solve a large DAE system. The DAE is obtained by > applying finite FV method to 3-D multi-specie ion/neutral fluids equations > with magnetic induction equation. The Jacobian for stiff part is formed by > using MatSetValuesStencil(). The Jacobian matrix is very sparse, no more than > 10 non-zeros on each row. MatSetValuesStencil requires local to global > mapping by calling ISLocaltoGlobalMapping(). Could you please instruct me how > to use local to global mapping? > > DMDA automatically sets up the ISLocaltoGlobalMapping() so you should not > need to. > > I also tried using DMCreateMatrix() to create the Jacobian. While local to > global mapping is not necessary, the matrix takes too much memory and > requires 64-bit indices. I would prefer to take the advantage of sparsity of > the Jacobian, pre-allocate the matrix to use as less as possible memory so > that the code can be run on a multi-core desktop. > > If using matrix-free does not work for you because the linear solves do not > converge or converge too slowly. Then you might be able to decrease the > memory used in the matrix. The DMCreateMatrix() does take advantage of > sparsity and tries to preallocate only what it needs. Often what it > preallocates is the best possible, but for multicomponent problems it has to > assume there is full coupling within all the degrees of freedom that live at > each at each grid point. How many components live at each grid point in your > model and are they not all coupled to each other in the equations? If they > are not fully coupled you can use either DMDASetBlockFills() or > DMDASetBlockFillsSparse() to indicate the reduced coupling that actually > exists for you model. > > Barry > > > > > > > > > Thank you very much for your advice. > > Jiannan