[petsc-users] Using the PCASM interface to define minimally overlapping subdomains

[Barry Smith]

On Sep 17, 2014, at 3:03 PM, Patrick Sanan <patrick.sanan at gmail.com> wrote:

> On 9/16/14 9:43 PM, Barry Smith wrote:
>> On Sep 16, 2014, at 2:29 PM, Matthew Knepley <knepley at gmail.com> wrote:
>> 
>>> On Tue, Sep 16, 2014 at 2:23 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>> 
>>>    Patrick,
>>> 
>>>      This "local part of the subdomains for this processor” term in  PCASMSetLocalSubdomains is, IMHO, extremely confusing. WTHWTS? Anyways, I think that if you set the is_local[] to be different than the is[] you will always end up with a nonsymetric preconditioner. I think for one dimension you need to use
>>> 
>>> No I don't think that is right. The problem below is that you have overlap in only one direction. Process 0 overlaps
>>> Process 1, but Process 1 has no overlap of Process 0. This is not how Schwarz is generally envisioned.
>>   Sure it is.
>>> Imagine the linear algebra viewpoint, which I think is cleaner here. You partition the matrix rows into non-overlapping
>>> sets. These sets are is_local[]. Then any information you get from another domain is another row, which is put into
>>> is[]. You can certainly have a non-symmetric overlap, which you have below, but it mean one way information
>>> transmission which is strange for convergence.
>>   No, not a all.
>> 
>> 
>> |    0      1      2      3     4      5     6    |
>> 
>>   Domain 0 is the region from  |   to  4  with Dirichlet boundary conditions at each end (| and 4). Domain 1 is from 2 to | with Dirichlet boundary conditions at each end (2 and |) .
>> 
>>   If you look at the PCSetUp_ASM() and PCApply_ASM() you’ll see all kinds of VecScatter creations from the various is and is_local, “restriction”, “prolongation” and “localization” then in the apply the different scatters are applied in the two directions, which results in a non-symmetric operator.
> 
> I was able to get my uniprocessor example to give the (symmetric) preconditioner I expected by commenting out the check in PCSetUp_ASM (line 311 in asm.c)

    if (firstRow != lastRow) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Specified ASM subdomain sizes were invalid: %d != %d", firstRow, lastRow);

   This check is absurd and needs to be removed. 

> and using PCASMSetLocalSubdomains with the same (overlapping) IS's for both is and is_local ([0 1 2 3] and [3 4 5 6] in the example above). It also works passing NULL for is_local.

  Great.
> 
> I assume that the purpose of the check mentioned above is to ensure that every grid point is assigned to exactly one processor, which is needed by whatever interprocess scattering goes on in the implementation. Also, I assume that augmenting the domain definition with an explicit specification of the way domains are distributed over processes allows for more controllable use of PC_ASM_RESTRICT, with all its attractive properties.
> 
> Anyhow, Barry's advice previously in this thread works locally (for one test case) if you remove the check above, but the current implementation enforces something related to what Matt describes, which might be overly restrictive if multiple domains share a process.  The impression I got initially from the documentation was that if one uses PC_ASM_BASIC, the choice of is_local should only influence the details of the communication pattern, not (in exact arithmetic, with process-count-independent subsolves) the preconditioner being defined.

   The “communication pattern” does determine the preconditioner being defined.

    The introduction of is_local[] broke the clean usage of PC_ASM_* that use to exist, so your confusion is our fault, not yours.
> 
> 
> For regular grids this all seems pretty pathological (in practice I imagine people want to use symmetric overlaps,

   As I said above you are using "symmetric overlaps,”. It just looks “unsymmetric” if you introduce this concept of “non-overlapping initial subdomains” which is an unneeded harmful concept.

   Barry


> and I assume that one domain per node is the most common use case), but I could imagine it being more of a real concern when working with unstructured grids.
> 
>> 
>>    Barry
>> 
>> 
>> 
>>>   Matt
>>>  
>>>> is[0] <-- 0 1 2 3
>>>> is[1] <-- 3 4 5 6
>>>> is_local[0] <-- 0 1 2 3
>>>> is_local[1] <-- 3 4 5 6
>>> Or you can pass NULL for is_local use PCASMSetOverlap(pc,0);
>>> 
>>>   Barry
>>> 
>>> 
>>> Note that is_local[] doesn’t have to be non-overlapping or anything.
>>> 
>>> 
>>> On Sep 16, 2014, at 10:48 AM, Patrick Sanan <patrick.sanan at gmail.com> wrote:
>>> 
>>>> For the purposes of reproducing an example from a paper, I'd like to use PCASM with subdomains which 'overlap minimally' (though this is probably never a good idea in practice).
>>>> 
>>>> In one dimension with 7 unknowns and 2 domains, this might look like
>>>> 
>>>> 0  1  2  3  4  5  6  (unknowns)
>>>> ------------          (first subdomain  : 0 .. 3)
>>>>         -----------  (second subdomain : 3 .. 6)
>>>> 
>>>> The subdomains share only a single grid point, which differs from the way PCASM is used in most of the examples.
>>>> 
>>>> In two dimensions, minimally overlapping rectangular subdomains would overlap one exactly one row or column of the grid. Thus, for example, if the grid unknowns were
>>>> 
>>>> 0  1  2  3  4  5  |
>>>> 6  7  8  9  10 11 | |
>>>> 12 13 14 15 16 17   |
>>>>         --------
>>>> -----------
>>>> 
>>>> then one minimally-overlapping set of 4 subdomains would be
>>>> 0 1 2 3 6 7 8 9
>>>> 3 4 5 9 10 11
>>>> 6 7 8 9 12 13 14 15
>>>> 9 10 11 15 16 17
>>>> as suggested by the dashes and pipes above. The subdomains only overlap by a single row or column of the grid.
>>>> 
>>>> My question is whether and how one can use the PCASM interface to work with these sorts of decompositions (It's fine for my purposes to use a single MPI process). In particular, I don't quite understand if should be possible to define these decompositions by correctly providing is and is_local arguments to PCASMSetLocalSubdomains.
>>>> 
>>>> I have gotten code to run defining the is_local entries to be subsets of the is entries which define a partition of the global degrees of freedom*, but I'm not certain that this was the correct choice, as it appears to produce an unsymmetric preconditioner for a symmetric system when I use direct subdomain solves and the 'basic' type for PCASM.
>>>> 
>>>> * For example, in the 1D example above this would correspond to
>>>> is[0] <-- 0 1 2 3
>>>> is[1] <-- 3 4 5 6
>>>> is_local[0] <-- 0 1 2
>>>> is_local[1] <-- 3 4 5 6
>>>> 
>>>> 
>>>> 
>>>> 
>>> 
>>> 
>>> 
>>> -- 
>>> 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