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