April 2011 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] ML changing the number of levels
by Leo van Kampenhout 20 Apr '11

20 Apr '11

Hello, I was wondering if it is already possible to set the number of levels with the ML package. I tried -pc_mg_levels and PCMGSetLevels, both do not work. I have an 8*10^6 Poisson matrix, and so far ML has given me the best results. However, I think with increasing the gridsize, i should also increase the number of multigrid levels. What are your thoughts? Thanks, Leo 2011/3/23 Barry Smith <bsmith(a)mcs.anl.gov> > > On Mar 22, 2011, at 9:42 PM, Jed Brown wrote: > > > I was hoping you would just fix it. :) > > I don't know how to fix it except with a hack. I'll look at as soon as > my brain grows back, hopefully in a day or two. > > > Barry > > > > > > >> On Mar 23, 2011 3:39 AM, "Barry Smith" <bsmith(a)mcs.anl.gov> wrote: > >> > >> > >> > >> On Mar 22, 2011, at 7:34 PM, Jed Brown wrote: > >> > >> > I just fixed a few things in PCReset_ML, but t... > >> > >> Right. I forgot about this issue. Maybe put a note about this in the > code for the future. > >> > >> > >> Barry > >> > > > >

1 0

Re: [petsc-dev] [petsc-users] FETI-DP
by Jed Brown 20 Apr '11

20 Apr '11

Thomas, we should move this discussion to petsc-dev, are you subscribed to that list? On Wed, Apr 20, 2011 at 13:55, Thomas Witkowski < thomas.witkowski(a)tu-dresden.de> wrote: > There one small thing on the implementation details of the FETI-DP, I > cannot figure out. Maybe some of you could help me to understand it, though > it is not directly related to PETSc. Non of the publications says something > about how to distribute the Lagrange multipliers over the processors. Is > there any good way to do it or can it done arbitrarily? > All their work that I have seen assumes a fully redundant set of Lagrange multipliers. In that context, each Lagrange multiplier only ever couples two subdomains together. Either process can then take ownership of that single Lagrange multiplier. > And should be the jump operators B^i be directly assembled or should they > be implemented in a matrix-free way? > Usually these constraints are sparse so I think it is no problem to assume that they are always assembled. > I'm confuse because in the work of Klawoon/Rheinbach, it is claimed that > the following operator can be solved in a pure local way: > > F = \sum_{i=1}^{N} B^i inv(K_BB^i) trans(B^i) > Did they use "F" for this thing? Usually F is the FETI-DP operator which involves a Schur complement of the entire partially assembled operator in the dual space. In any case, this thing is not purely local since the jump operators B^i need neighboring values so it has the same communication as a MatMult. > With B^i the jump operators and K_BB^i the discretization of the sub > domains with the primal nodes. > I think you mean "with the primal nodes removed". > From the notation it follows that EACH local solve takes the whole vector > of Lagrange multipliers. But this is not applicable for a good parallel > implementation. Any hint on this topic would be helpful for me to understand > this problem. > I can't tell from their papers how B is stored. It would be natural to simply store B as a normal assembled matrix with a standard row partition of the Lagrange multipliers. Then you would apply the subdomain solve operator using MatMultTranspose(B,XLambdaGlobal,XGlobal); for (i=0; i<nlocalsub; i++) { Vec XSubdomain,YSubdomain; VecGetSubVector(XGlobal,sublocal[i],&XSubdomain); // no copy if subdomains are contiguous VecGetSubVector(YGlobal,sublocal[i],&YSubdomain); // also no copy KSPSolve(kspK_BB[i],XSubdomain,YSubdomain); // purely local solve, often KSPPREONLY and PCLU VecRestoreSubVector(XGlobal,sublocal[i],&XSubdomain); VecRestoreSubVector(YGlobal,sublocal[i],&YSubdomain); } MatMult(B,YGlobal,YLambdaGlobal); All the communication is in the MatMultTranspose and MatMult. The "Global" vectors here are global with respect to K_BB (interior and interface dofs, primal dofs removed). I don't think there is need to ever store K_BB as a parallel matrix, it would be a separate matrix per subdomain (in the general case, subdomains could be parallel on subcommunicators). This code should handle nlocalsub subdomains owned by the local communicator, typically PETSC_COMM_SELF. The index sets (IS) in sublocal represent the global (space of K_BB) dofs, usually these are contiguous sets so they can be represented very cheaply. Barry, would you do it differently?

1 0

Re: [petsc-dev] DACreate3d should work for 1 node in Z
by Matthew Knepley 19 Apr '11

19 Apr '11

On Tue, Apr 19, 2011 at 6:19 PM, Jed Brown <jed(a)59a2.org> wrote: > On Wed, Apr 20, 2011 at 00:35, Ethan Coon <ecoon(a)lanl.gov> wrote: > >> Ok, here's a patch that does the very very special case of: >> >> p=P=1, s > p, DMDA_BOUNDARY_PERIODIC in the z-direction. >> >> Note this is really only possible in the z-direction... >> > > That is the case that you usually want, otherwise DMDAVecGetArray has > rather inefficient indexing. > > >> when in the >> z-direction, we can just check if the global index < 0 or > x*y*z and >> adjust appropriately. So this can't generalize to x- or y- directions >> (though it could do the y-direction in the 2D case, allowing one to do >> 1D problems in a 2D algorithm?) >> >> Compared to the rest of the cruft in da3.c and how non-general the >> global index generation is overall, it's not actually that ugly... >> > > We should all go unstructured for everything. Then we would have already > paid the storage cost of heavier data structures so many of the special > cases go away. :-) > > I did this already in src/dm/impls/cartesian, but we don't want to depend on Boost for the data structure. I will rewrite this with a simple multiindex container in PETSc when I get a chance. Once this is done, the code is very very simple to construct the scatter. Matt > >> I've tested with s=2 and 3, BOX and STAR, and the indices look right >> (i.e. they are identical in the z-dimension). But it's ugly, so please >> test with your stuff too Jed. >> > > It seems to be working with Andrew's MHD code so I've pushed your > patch. Thank you Ethan. > > Andrew, your code is running now with NZ=1. Your function evaluation still > evaluates to zero so SNESSolve is not doing any iterations. I assume there > is some other step to turn on your source terms. > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

1 0

Re: [petsc-dev] DACreate3d should work for 1 node in Z
by Jed Brown 19 Apr '11

19 Apr '11

On Wed, Apr 20, 2011 at 00:35, Ethan Coon <ecoon(a)lanl.gov> wrote: > Ok, here's a patch that does the very very special case of: > > p=P=1, s > p, DMDA_BOUNDARY_PERIODIC in the z-direction. > > Note this is really only possible in the z-direction... > That is the case that you usually want, otherwise DMDAVecGetArray has rather inefficient indexing. > when in the > z-direction, we can just check if the global index < 0 or > x*y*z and > adjust appropriately. So this can't generalize to x- or y- directions > (though it could do the y-direction in the 2D case, allowing one to do > 1D problems in a 2D algorithm?) > > Compared to the rest of the cruft in da3.c and how non-general the > global index generation is overall, it's not actually that ugly... > We should all go unstructured for everything. Then we would have already paid the storage cost of heavier data structures so many of the special cases go away. :-) > > I've tested with s=2 and 3, BOX and STAR, and the indices look right > (i.e. they are identical in the z-dimension). But it's ugly, so please > test with your stuff too Jed. > It seems to be working with Andrew's MHD code so I've pushed your patch. Thank you Ethan. Andrew, your code is running now with NZ=1. Your function evaluation still evaluates to zero so SNESSolve is not doing any iterations. I assume there is some other step to turn on your source terms.

1 0

Re: [petsc-dev] DACreate3d should work for 1 node in Z
by Ethan Coon 19 Apr '11

19 Apr '11

Ok, here's a patch that does the very very special case of: p=P=1, s > p, DMDA_BOUNDARY_PERIODIC in the z-direction. Note this is really only possible in the z-direction... when in the z-direction, we can just check if the global index < 0 or > x*y*z and adjust appropriately. So this can't generalize to x- or y- directions (though it could do the y-direction in the 2D case, allowing one to do 1D problems in a 2D algorithm?) Compared to the rest of the cruft in da3.c and how non-general the global index generation is overall, it's not actually that ugly... I've tested with s=2 and 3, BOX and STAR, and the indices look right (i.e. they are identical in the z-dimension). But it's ugly, so please test with your stuff too Jed. Ethan On Tue, 2011-04-19 at 23:28 +0200, Jed Brown wrote: > On Tue, Apr 19, 2011 at 23:21, Ethan Coon <ecoon(a)lanl.gov> wrote: > The attached patch fixes the errors to error in the correct > cases. > > Pushed, thanks. -- ------------------------------------ Ethan Coon Post-Doctoral Researcher Applied Mathematics - T-5 Los Alamos National Laboratory 505-665-8289 http://www.ldeo.columbia.edu/~ecoon/ ------------------------------------

1 0

Re: [petsc-dev] DACreate3d should work for 1 node in Z
by Jed Brown 19 Apr '11

19 Apr '11

On Tue, Apr 19, 2011 at 23:21, Ethan Coon <ecoon(a)lanl.gov> wrote: > The attached patch fixes the errors to error in the correct cases. Pushed, thanks.

1 0

Re: [petsc-dev] DACreate3d should work for 1 node in Z
by Ethan Coon 19 Apr '11

19 Apr '11

The attached patch fixes the errors to error in the correct cases. Ethan On Tue, 2011-04-19 at 22:48 +0200, Jed Brown wrote: > On Tue, Apr 19, 2011 at 22:44, Ethan Coon <ecoon(a)lanl.gov> wrote: > The short answer -- make the z-size = s. It's not that huge > of a hit -- > the local arrays are 6*NX*NY*ndofs instead of 5*NX*NY*ndofs. > > > > Indeed, though it is potentially about twice as much work because you > only have to evaluate residuals on the owned nodes. > > > Otherwise, I'll try to implement the special-special case of > DMDA_BOUNDARY_PERIODIC, P=1, zl=1, and s > 1 soon. > > Great, thanks. -- ------------------------------------ Ethan Coon Post-Doctoral Researcher Applied Mathematics - T-5 Los Alamos National Laboratory 505-665-8289 http://www.ldeo.columbia.edu/~ecoon/ ------------------------------------

1 0

Re: [petsc-dev] DACreate3d should work for 1 node in Z
by Jed Brown 19 Apr '11

19 Apr '11

On Tue, Apr 19, 2011 at 22:44, Ethan Coon <ecoon(a)lanl.gov> wrote: > The short answer -- make the z-size = s. It's not that huge of a hit -- > the local arrays are 6*NX*NY*ndofs instead of 5*NX*NY*ndofs. > Indeed, though it is potentially about twice as much work because you only have to evaluate residuals on the owned nodes. > > Otherwise, I'll try to implement the special-special case of > DMDA_BOUNDARY_PERIODIC, P=1, zl=1, and s > 1 soon. > Great, thanks.

1 0

Re: [petsc-dev] DACreate3d should work for 1 node in Z
by Ethan Coon 19 Apr '11

19 Apr '11

On Tue, 2011-04-19 at 22:25 +0200, Jed Brown wrote: > On Tue, Apr 19, 2011 at 22:19, Ethan Coon <ecoon(a)lanl.gov> wrote: > At any rate, I'll fix the error message in 2D and add it to 3D > to only > error in the cases that actually currently break. Doing > DMDA_BOUNDARY_PERIODIC with P=1 and s > zl (the multiple-wraps > around > the domain case) will take more thought... > > The special case of 1 is probably the most important since it is the > minimal size for a 3D method running a 2D problem. Agreed -- I was surprised this didn't work because we've been doing a problem where we use stencil_size > 1, and we run z-dimension size = 1 problems all the time for this case. Apparently we only do stencil_size > 1 in the 2D case (but intended to generalize to 3D soon). The periodic case is the one you need to collapse 3D into 2D with no extra work... an infinite domain in the 3rd dimension. The short answer -- make the z-size = s. It's not that huge of a hit -- the local arrays are 6*NX*NY*ndofs instead of 5*NX*NY*ndofs. Otherwise, I'll try to implement the special-special case of DMDA_BOUNDARY_PERIODIC, P=1, zl=1, and s > 1 soon. Ethan > Multiple wraps gets hairy pretty quickly and it doesn't seem important > to support in the short term. It's certainly less important than > general 3D domains with some bits empty (since this comes up > frequently for coarse levels of multigrid). -- ------------------------------------ Ethan Coon Post-Doctoral Researcher Applied Mathematics - T-5 Los Alamos National Laboratory 505-665-8289 http://www.ldeo.columbia.edu/~ecoon/ ------------------------------------

1 0

Re: [petsc-dev] DACreate3d should work for 1 node in Z
by Jed Brown 19 Apr '11

19 Apr '11

On Tue, Apr 19, 2011 at 22:19, Ethan Coon <ecoon(a)lanl.gov> wrote: > At any rate, I'll fix the error message in 2D and add it to 3D to only > error in the cases that actually currently break. Doing > DMDA_BOUNDARY_PERIODIC with P=1 and s > zl (the multiple-wraps around > the domain case) will take more thought... > The special case of 1 is probably the most important since it is the minimal size for a 3D method running a 2D problem. Multiple wraps gets hairy pretty quickly and it doesn't seem important to support in the short term. It's certainly less important than general 3D domains with some bits empty (since this comes up frequently for coarse levels of multigrid).

1 0