On Sat, Jan 15, 2011 at 8:23 AM, Gianluca Meneghello <gianmail at gmail.com>wrote:
> Dear Barry > > In order to use matrix-free methods with PCFIELDSPLIT, is it correct > to start from > > MatShellSetOperation(mat,MATOP_MATGETSUBMATRIX, (void(*)(void)) > PetscErrorCode (*UserMatGetSubMatrix)(......)); > > or does it require to dive deeper inside PETSc code? > That is fine. You must return a Mat object, probably a MatShell in your case, which applys the requested block of the operator. Matt > Thanks > > Gianluca > > On 13 January 2011 22:16, Barry Smith <bsmith at mcs.anl.gov> wrote: > > > > On Jan 13, 2011, at 2:09 PM, Gianluca Meneghello wrote: > > > >> Barry, > >> > >> I'm afraid I didn't explain myself well in my first question. > >> > >> If I use PCFIELDSPLIT with > PCFieldSplitSetType(pc,PC_COMPOSITE_ADDITIVE), is it the same as using > PCASM? > > > > Same in what sense? It solves a bunch of subproblems independently and > adds together all the solutions. There can be overlapping in the fields or > not depending how what you choose. The decomposition in ASM is by "geometry" > while the decomposition in the PCFIELDSPLIT is between different "fields" or > "types of variables". So yes they have many similarities. > > > > > >> > >> Concerning the MATSHELL to be used with PCFIELDSPLIT, is there an > example from where to start from? I guess I'm one of the people for which > it is non-trivial :-) > > > > Not really. > > > > Barry > > > >> > >> Thanks > >> > >> Gianluca > >> > >> Il giorno 13/gen/2011, alle ore 20.14, Barry Smith ha scritto: > >> > >>> > >>> On Jan 13, 2011, at 12:00 PM, Gianluca Meneghello wrote: > >>> > >>>> Dear Barry and Jed, > >>>> > >>>> thanks again for your answers. I'm at the moment trying to understand > more about how ASM and FieldSplit works. I've started reading Barry's book > "Domain Decomposition, Parallel Multilevel Methods for Elliptic Partial > Differential Equations". I hope it's the right starting point. > >>>> > >>>> Please let me know if you have further suggested readings, taking into > account that I know nothing on domain decomposition in general --- but > something in multigrid. > >>>> > >>>> Let me ask you a couple of questions for the moment: > >>>> Is PCFIELDSPLIT additive the same as PCASM (I guess no as it does not > use overlap)? > >>> > >>> It can be additive or multiplicative depending on what you set with > PCFieldSplitSetType() > >>> > >>>> Or is it a Substructuring Method (I haven't yet arrived to that > chapter of the book!). > >>> > >>> No, nothing to do with substructuring. > >>> > >>>> What's the difference between multiplicative and > symmetric-multiplicative for PCFIELDSPLIT? Does symmetric-multiplicative > refers to eq 1.15 of "Domain Decomposition"? > >>> > >>> Yes > >>> > >>>> > >>>> Is it possible use ASM and/or FieldSplit with a matrix-free method? > >>> > >>> Yes, BUT the algorithms are coded around MatGetSubMatrix() and or > MatGetSubMatrices() so to do matrix free you need to have code that applies > "part" of the operator at a time (that is you cannot just have a matrix > vector product that applies the entire operator to the entire vector. Once > you have the ability to apply "part" of the operator at a time you need to > code up a MATSHELL that responds appropriately to MatGetSubMatrix() and or > MatGetSubMatrices() and returns new matrix-free shell matrices that apply > only "their" part of the operator. This is non-trivial for many people but > possible. > >>> > >>> Barry > >>> > >>> > >>>> > >>>> Thanks > >>>> > >>>> Gianluca > >>>> > >>>> Il giorno 06/gen/2011, alle ore 21.55, Barry Smith ha scritto: > >>>> > >>>>> > >>>>> On Jan 6, 2011, at 6:01 AM, Gianluca Meneghello wrote: > >>>>> > >>>>>> Dear Barry, > >>>>>> > >>>>>> thanks a lot for your answer. > >>>>>> > >>>>>> I tried to do some experiments with MatGetSubMatrix, but I guess I'm > >>>>>> doing something wrong as it looks like being painfully slow. > >>>>> > >>>>> Hmm, we've always found the getsubmatrix takes a few percent of the > time. Perhaps you are calling it repeatedly for the each domain, rather than > once and reusing it? Also use MatGetSubMatrices() and get all the > submatrices in one call rather than one at a time. > >>>>> > >>>>>> > >>>>>> I changed approach and now I'm using the ASM preconditioner. What > I'm > >>>>>> actually trying to do is to split the domain in different parts --- > >>>>>> like interior and boundaries --- and relax (solve) each one with a > >>>>>> different smoother (solver). In your opinion, is this the right > >>>>>> approach? > >>>>> > >>>>> Worth trying since it is easy. You can experiment with different > smoothers on the subdomains using the -sub_pc_type etc options and set > different prefixes for different subdomains. > >>>>> > >>>>>> So far it looks much faster than my previous approach of > >>>>>> extracting each submatrix. > >>>>> > >>>>> ASM just uses MatGetSubMatrices() so shouldn't be faster or slower > than a custom code that does the same thing. > >>>>> > >>>>>> > >>>>>> Also, please let me ask you one more thing. When using ASM with > >>>>>> different subdomains on the same process, is the order in which the > >>>>>> domains are solved the same as the one in which they are stored in > the > >>>>>> IS array passed to PCASMSetLocalSubdomains()? I would be interested > in > >>>>>> controlling this in order to build a downstream marching smoother. > >>>>> > >>>>> It is only additive, there is no order as Jed noted. Doing > multiplicative in general is tricky because you want to just update the > parts of the residual that need to be updated. > >>>>> > >>>>>> > >>>>>> Looking at the references, I've noticed you have worked on > multigrid. > >>>>>> What I'm trying to do is close to what is described in Diskin, > Thomas, > >>>>>> Mineck, "Textbook Multigrid Efficiency for Leading Edge Stagnation", > >>>>>> in case you already know the paper. > >>>>>> > http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20040081104_2004082284.pdf > >>>>>> > >>>>>> Again, thanks a lot. > >>>>>> > >>>>>> Gianluca > >>>>>> > >>>>>> On 3 January 2011 17:43, Barry Smith <bsmith at mcs.anl.gov> wrote: > >>>>>>> > >>>>>>> Gianluca, > >>>>>>> > >>>>>>> The expected use is with the VecScatter object. First you create a > VecScatter object with VecScatterCreate() then each time you need the > "subvector" you call VecScatterBegin() followed by VecScatterEnd() Note that > usually the VecScatter object is retained and used many times. > >>>>>>> > >>>>>>> Barry > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>>> On Jan 3, 2011, at 5:22 AM, Gianluca Meneghello wrote: > >>>>>>> > >>>>>>>> Hi, > >>>>>>>> > >>>>>>>> I'm new to PETSc, so that this can be a very simple question: > >>>>>>>> > >>>>>>>> I'm looking for something like VecGetSubVector, which I've seen it > >>>>>>>> exists in the dev version but not in the released one. > >>>>>>>> > >>>>>>>> I need to write a smoother for a multigrid algorithm (something > like a > >>>>>>>> block Gauss Seidel) which can be written in matlab as > >>>>>>>> > >>>>>>>> for j = 1:ny > >>>>>>>> P = <some matrix indices as function of j>; > >>>>>>>> du(P) = L(P,P) \ ( rhs(P) - L(P,:)*du + D2(P,P)*du(P) ); > >>>>>>>> end > >>>>>>>> > >>>>>>>> where L is a matrix (in my case the linearized Navier Stokes). > >>>>>>>> > >>>>>>>> I was thinking about using IS for declaring P, so that D2(P,P) can > be > >>>>>>>> obtained using MatGetSubMatrix. I would need the same for the > vector > >>>>>>>> du. > >>>>>>>> > >>>>>>>> Is there a way to do that without using the developer version? (I > >>>>>>>> really don't feel like being "experienced with building, using and > >>>>>>>> debugging PETSc). > >>>>>>>> > >>>>>>>> Thanks in advance > >>>>>>>> > >>>>>>>> Gianluca > >>>>>>> > >>>>>>> > >>>>> > >>>> > >>> > >> > > > > > -- 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 -------------- next part -------------- An HTML attachment was scrubbed... URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110115/b809ebcb/attachment-0001.htm>
