[petsc-dev] DMPlex + P_1 FE interpolation

Matthew Knepley knepley at gmail.com
Mon Feb 10 07:35:09 CST 2020 

On Mon, Feb 10, 2020 at 5:31 AM Dave May <dave.mayhem23 at gmail.com> wrote:

>
>
> On Mon 10. Feb 2020 at 14:17, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Mon, Feb 10, 2020 at 5:08 AM Dave May <dave.mayhem23 at gmail.com> wrote:
>>
>>>
>>>
>>> On Mon 10. Feb 2020 at 13:09, Matthew Knepley <knepley at gmail.com> wrote:
>>>
>>>> On Sun, Feb 9, 2020 at 11:11 PM Pierre Jolivet <
>>>> pierre.jolivet at enseeiht.fr> wrote:
>>>>
>>>>>
>>>>>
>>>>> On 10 Feb 2020, at 6:20 AM, Matthew Knepley <knepley at gmail.com> wrote:
>>>>>
>>>>> On Sun, Feb 9, 2020 at 3:23 PM Pierre Jolivet <
>>>>> pierre.jolivet at enseeiht.fr> wrote:
>>>>>
>>>>>> Hello,
>>>>>> I’ve a hard time answering the following DMPlex questions by just
>>>>>> looking at some of the examples and manual.
>>>>>> Considering two DMPlex dm and dma, as in
>>>>>> petsc/src/dm/impls/plex/examples/tests/ex19.c, I’d like to interpolate a
>>>>>> simple P_1 FE function from dm to dma.
>>>>>> The DMCreateInterpolation call gives me:
>>>>>> [0]PETSC ERROR: Invalid argument
>>>>>> [0]PETSC ERROR: Number of fine indices 0 != 4 dual basis vecs
>>>>>> […]
>>>>>>
>>>>>
>>>>> It looks like your fine grid has no discretization, since 0 is
>>>>> numFIndices from
>>>>>
>>>>>   ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell,
>>>>> &numFIndices, &findices, NULL);CHKERRQ(ierr);
>>>>>
>>>>>
>>>>>> [0]PETSC ERROR: #1 DMPlexComputeInterpolatorGeneral() line 2508 in
>>>>>> petsc/src/dm/impls/plex/plexfem.c
>>>>>> [0]PETSC ERROR: #2 DMCreateInterpolation_Plex() line 7688 in
>>>>>> petsc/src/dm/impls/plex/plex.c
>>>>>> [0]PETSC ERROR: #3 DMCreateInterpolation() line 1139 in
>>>>>> petsc/src/dm/interface/dm.c
>>>>>> But the DMs look OK, don’t they, cf. below?
>>>>>> So I have three simple questions:
>>>>>> 1) are all tests at the bottom of ex19.c broken because of PRAgMaTIc
>>>>>> or because of DMPlex currently not supporting some operations? (I’m not
>>>>>> using PRAgMaTIc to do mesh adaptation, so I was hoping to not run into an
>>>>>> error)
>>>>>>
>>>>>
>>>>> I don't think its broken.
>>>>>
>>>>>
>>>>> Oh, OK. Could you help me figure out what’s the problem then, e.g.,
>>>>> with a slight (command line) variation of test #6, please?
>>>>>
>>>>
>>>> Sure. I am at SIAM this week, but as soon as I can I will get you the
>>>> fix.
>>>>
>>>>
>>>>> $ cd src/dm/impls/plex/examples/tests
>>>>> $ git diff ex19.c
>>>>> $ make ex19
>>>>> $ mpirun ./ex19 -dim 3 -nbrVerEdge 10 -dm_plex_separate_marker 0 -met
>>>>> 0 -hmin 0.1 -hmax 0.3 -init_dm_view -adapt_dm_view -do_L2
>>>>> -petscspace_degree 1 -petscfe_default_quadrature_order 1
>>>>> -dm_plex_hash_location
>>>>> [0]PETSC ERROR: Nonconforming object sizes
>>>>> [0]PETSC ERROR: The section point closure size 0 != dual space
>>>>> dimension 4
>>>>> […]
>>>>> [0]PETSC ERROR: #1 DMProjectLocal_Generic_Plex() line 633 in
>>>>> src/dm/impls/plex/plexproject.c
>>>>> [0]PETSC ERROR: #2 DMProjectFunctionLocal_Plex() line 771 in
>>>>> src/dm/impls/plex/plexproject.c
>>>>> [0]PETSC ERROR: #3 DMProjectFunctionLocal() line 7809 in
>>>>> src/dm/interface/dm.c
>>>>> [0]PETSC ERROR: #4 DMProjectFunction() line 7766 in
>>>>> src/dm/interface/dm.c
>>>>>
>>>>> If I comment the DMProjectFunction() call, I end up with the same
>>>>> error as in my first message in DMCreateInterpolation().
>>>>>
>>>>> 2) is DMCreateInterpolation + MatInterpolate the correct way of
>>>>>> transferring one Vec from a DMPlex onto another?
>>>>>>
>>>>>
>>>>> That is the intent.
>>>>>
>>>>>
>>>>>> 3) if yes, by looking at the names of the arguments in
>>>>>> DMPlexComputeInterpolatorGeneral, dmc and dmf, could you comment on the
>>>>>> performance of this function for unrelated meshes, e.g., if both DMs are
>>>>>> “fine” and not one coarse and the other fine (albeit non-nested), for
>>>>>> simple P_k spaces.
>>>>>>
>>>>>
>>>>> In general, it is going to be horrible. Here is what it does: locate
>>>>> the fine quadrature points in the coarse grid and interpolate to them. This
>>>>> quadrature can have huge errors if it falls across multiple cells. This is
>>>>> why the nested version works perfectly, and also why Patrick Farrell and
>>>>> James Maddison have the Supermesh library, which makes a refinement of the
>>>>> mesh until the quadrature is accurate everywhere. That way they guarantee
>>>>> that at least the zeroth moment is preserved.
>>>>>
>>>>>
>>>>> Two subquestions if I may:
>>>>> 1) are there any plans to have this integrated through an external
>>>>> package?
>>>>>
>>>>
>>>> In the absence of a plan, there is a hope. I would really like it to
>>>> happen.
>>>>
>>>>
>>>>> 2) if I understand you correctly, you answered about the numerical
>>>>> performance of the function. I can live with high interpolation errors if
>>>>> both meshes are “far" from each other. I was mostly interested in the
>>>>> parallel performance of the function.
>>>>>
>>>>
>>>> Everything is purely local except for point location. Since it has
>>>> never really been tested in this mode, I am sure the
>>>> scaling can be bad. I believe the default is to extrapolate if the
>>>> point is not covered, which makes sense for mostly
>>>> matching meshes. There is parallel point location, but it is intended
>>>> for a few points where we are sampling the solution,
>>>> rather than lots of points everywhere which you would get for
>>>> non-matching meshes with different distributions. Could
>>>> you say what kind of situation you are trying to optimize for?
>>>>
>>>
>>> Matt, probably clarifying the parallel point location algorithm is
>>> helpful.
>>>
>>> 1/ Does it broadcast all off rank points to every rank?
>>> Or
>>> 2/ Does it broadcast sub domain bounding boxes  to every rank, and then
>>> scatter points to candidate owning ranks based on the bounding boxes
>>> containing off rank points?
>>>
>>> I recall the method (years ago) did  what’s described in (1)
>>>
>>
>> I thought it was 2.
>>
>> This is one of the things on the list when we scale up the plasmas
>> physics code. My reading
>> says that a hierarchy of regular divisions is just as efficient for
>> location as a structure for
>> irregular divisons, like a k-d tree, and we can keep the bottom part of
>> the hierarchy everywhere.
>> Is that your feeling?
>>
>
> Yes.
>
>  (2) is good for general purpose usage. I believe anything better requires
> restrictive assumptions which render the interpolation suitable only for
> specific use case (eg assuming  an off rank point lives within a nearest or
> next nearest neighbor sub domain)
>

We have to work a little harder, since with a hierarchy you are sending
candidates to a subgroup of processors for a more local
search, so you probably don't want to send everything to some "master"
process, but split it up among them. This is exactly the
message aggregation problem again.

   Matt


>
>>   Matt
>>
>>
>>>
>>>>   Thanks,
>>>>
>>>>     Matt
>>>>
>>>>
>>>>> Thanks,
>>>>> Pierre
>>>>>
>>>>>   Thanks,
>>>>>
>>>>>     Matt
>>>>>
>>>>>
>>>>>> Thanks in advance for your help,
>>>>>> Pierre
>>>>>>
>>>>>> $ mpirun -n 1 ./ex19 -msh in.msh -init_dm_view ::ascii_info
>>>>>> -adapt_dm_view ::ascii_info -mat_view ::ascii_info -do_L2
>>>>>> -petscspace_degree 1
>>>>>> DM Object: DMinit 1 MPI processes
>>>>>>   type: plex
>>>>>> DMinit in 3 dimensions:
>>>>>>   0-cells: 1331
>>>>>>   1-cells: 7930
>>>>>>   2-cells: 12600
>>>>>>   3-cells: 6000
>>>>>> Labels:
>>>>>>   depth: 4 strata with value/size (0 (1331), 1 (7930), 2 (12600), 3
>>>>>> (6000))
>>>>>>   Face Sets: 6 strata with value/size (4 (200), 1 (200), 5 (200), 2
>>>>>> (200), 3 (200), 6 (200))
>>>>>>   Cell Sets: 1 strata with value/size (0 (6000))
>>>>>> DM Object: DMadapt (adapt_) 1 MPI processes
>>>>>>   type: plex
>>>>>> DMadapt in 3 dimensions:
>>>>>>   0-cells: 2905
>>>>>>   1-cells: 18888
>>>>>>   2-cells: 31368
>>>>>>   3-cells: 15384
>>>>>> Labels:
>>>>>>   depth: 4 strata with value/size (0 (2905), 1 (18888), 2 (31368), 3
>>>>>> (15384))
>>>>>>   Face Sets: 6 strata with value/size (1 (200), 4 (200), 6 (200), 2
>>>>>> (200), 5 (200), 3 (200))
>>>>>>   Cell Sets: 1 strata with value/size (0 (15384))
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>> https://www.cse.buffalo.edu/~knepley/
>>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>>
>>>>>
>>>>>
>>>>
>>>> --
>>>> 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
>>>>
>>>> https://www.cse.buffalo.edu/~knepley/
>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>
>>>
>>
>> --
>> 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
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>

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

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20200210/e1491997/attachment-0001.html>

More information about the petsc-dev mailing list