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? 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 >>>>>>> >>>>>>> >>>>> >>>> >>> >> > >
