[petsc-users] DMPlex Interpolation

Mon Aug 29 08:55:02 CDT 2022

On Sun, Aug 28, 2022 at 7:17 PM Matthew Knepley <knepley at gmail.com> wrote:

> On Sun, Aug 28, 2022 at 5:36 PM Mike Michell <mi.mike1021 at gmail.com>
> wrote:
>
>> Thank you for the reply.
>>
>> I think it can be more helpful for me if the attached sample code
>> (DMInterpolation_Mod.tar) could be checked by you.
>>
>
> Okay, you are right. I will run it tomorrow.
>

The problem was somewhat conceptual. When you are restricting state,
meaning a function you want to use directly, rather than
residuals, meaning some average or integral, you want to preserve scale.
This is what the 'vscale' vector is for. I added the
VecPointwiseMult to your code (attached) to do this. I also changed to a
more general mehs creation. To get your run, use

  ./App.x -dm_plex_filename square.msh

There is inaccuracy, but it is order O(h) just like the interpolation error
from constants. You can see this by refining the mesh

  ./App.x -dm_plex_simplex 0 -dm_plex_box_faces 25,25

I am attaching the picture. The last value is 0.98 instead of 1, but the
error is decreasing like h.

  Thanks,

     Matt

>   Thanks,
>
>      Matt
>
>
>> If you run the sample code, basically, you would get
>> "VertexField_Map_from_Cell.vtu" file which displays a solution vector field
>> defined on each vertex, which is mapped from another solution field defined
>> on each cell-centroid. The vertex solution field in the vtu file should
>> match with "x-coordinate" value in the same file. Basically, the attached
>> code tries to reproduce vertex's x-coordinate value using
>> cell-centroid's x-coordinate value by means of finite-element space mapping
>> between Q0 and Q1.
>> From the vtu file, it is obvious that at the physical boundaries (left,
>> right, top, and bottom) x-coordinate value is not properly mapped from
>> cell-centroid, whereas interior vertex's solution seems okay. What I do not
>> understand is why the boundary vertices cannot have properly mapped
>> coordinate values? I believe I am doing something wrong and I was
>> suspicious about the ghost cell, which I do not have. But your comment
>> seems to say that the problem is not related to the ghost cell. What am I
>> doing wrong?
>>
>> Does your section mark these unknowns as constrained?
>> --> Yes, as shown in the code, I defined one section for each dmplex, one
>> for vertex and one for cell, respectively, with 1-dof for each (See code
>> line from 91 to 124).
>>
>> You should not get 0 unless you constrained these variables in the
>> section. If not, you should see it use the constant values in the cell.
>> --> Yes, you are correct. Now I do not have zero value at the physical
>> boundaries, however, it is still not proper value and seems underscaled,
>> compared to the proper answer (i.e., compared with the properly mapped
>> interior values).
>>
>> Thanks,
>>
>>
>>> On Sun, Aug 28, 2022 at 12:11 PM Mike Michell <mi.mike1021 at gmail.com>
>>> wrote:
>>>
>>>> Thank you for the reply.
>>>>
>>>> *I cannot quite understand. Are you saying that you have inhomogeneous
>>>> Dirichlet conditions on the boundary, and a 0 guess in the interior, and
>>>> you get all zeros from interpolation? Yes, we have no way of getting
>>>> boundary values into the operator. You could conceivably make an operator
>>>> that mapped between local vectors, but it is hard to see why you would want
>>>> this. In the solver, we usually  just care about the unknowns we are
>>>> solving for. Did I understand your question, or is it about something else?*
>>>> --> It is not imposing a physical boundary condition (Dirichlet
>>>> condition for example), but my solver needs to import an external solution
>>>> vector field resolved to each cell-centroid, then I need to map it to a
>>>> vertex field in which my solver is resolved for solving governing
>>>> equations. Then, my code solves its equations. After my solver is done
>>>> solving the equations, the vertex field should be re-map into cell-field
>>>> and then push back the cell field to the external code. With this context,
>>>> I have zero values returned at the vertices on the physical boundaries
>>>>
>>>
>>> Does your section mark these unknowns as constrained?
>>>
>>>
>>>> if I try mapping from cell-center to vertex field without having ghost
>>>> cell layers inside physical boundaries. For interior vertices, I get
>>>> reasonable mapped values, although the mapping accuracy seems to be able to
>>>> be improved as I mentioned in the earlier email.
>>>>
>>>
>>> As I replied, I think this is not true.
>>>
>>>
>>>> If dmplex for cell-center field would be defined to have ghost cell
>>>> layer, and fill it (it meant import ghost cell values as well from the
>>>> external code), and map it to vertex-field, can it be possible to fill
>>>> vertex points on the physical boundaries?
>>>>
>>>
>>> You should not get 0 unless you constrained these variables in the
>>> section. If not, you should see it use the constant values in the cell.
>>>
>>>   Thanks,
>>>
>>>      Matt
>>>
>>>
>>>> I was trying to test this strategy, but it
>>>> seems DMPlexComputeCellGeometryFVM() gives NaN values if dmplex has ghost
>>>> cell by calling DMPlexConstructGhostCells(). In detail,
>>>> DMPlexComputeCellGeometryFVM() gives zero area for each cell, and returns
>>>> NaN values for cell centroids and normals for the ghost cells.
>>>>
>>>> Thanks,
>>>>
>>>>
>>>>> On Sun, Aug 28, 2022 at 8:51 AM Mike Michell <mi.mike1021 at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Hi, thank you for the reply.
>>>>>>
>>>>>> I was able to manage mapping from cell-center to vertex. Basically,
>>>>>> in Fortran, it seems DMCreateInterpolation() requires the optional scaling
>>>>>> vector as a mandatory argument, which is strange.
>>>>>>
>>>>>
>>>>> It needs a custom Fortran wrapper to check for the NULL input from
>>>>> Fortran.
>>>>>
>>>>>
>>>>>> My dmplex has zero overlap layer over the procs, and there is no
>>>>>> ghost cell inside physical boundary. In this case, it seems mapping between
>>>>>> cell center to vertex returns zero (in case the global cell vector
>>>>>> initialized to zero) value at the nodes, which are located at physical
>>>>>> boundary. To manage this problem, is it mandatory to have ghost cells. Is
>>>>>> this correct understanding?
>>>>>>
>>>>>
>>>>> I cannot quite understand. Are you saying that you have inhomogeneous
>>>>> Dirichlet conditions on the boundary, and a 0 guess in the interior, and
>>>>> you get all zeros from interpolation? Yes, we have no way of getting
>>>>> boundary values into the operator. You could conceivably make an operator
>>>>> that mapped between local vectors, but it is hard to see why you would want
>>>>> this. In the solver, we usually  just care about the unknowns we are
>>>>> solving for.
>>>>>
>>>>> Did I understand your question, or is it about something else?
>>>>>
>>>>>
>>>>>> Also, my mapping accuracy itself seems to be improved. For simple
>>>>>> test, cell-center's x-coordinate value of each cell mapped to node and
>>>>>> printed to the vertex field as .vtu format, and the mapped
>>>>>> x-vertex-coordinate is quite different with the actual nodal coordinate
>>>>>> values what PETSc intrinsically provides through DMGetCoordinatesLocal(). I
>>>>>> believe I am doing something wrong, probably arguments
>>>>>> in PetscFECreateLagrange() can improve the mapping accuracy in
>>>>>> finite-element space?
>>>>>>
>>>>>
>>>>> Yes, a cellwise field is much less accurate than a linear field.
>>>>> Pushing the values up to a linear field does not make them more accurate.
>>>>>
>>>>>   Thanks,
>>>>>
>>>>>      Matt
>>>>>
>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>>
>>>>>>> On Thu, Aug 25, 2022 at 7:12 PM Mike Michell <mi.mike1021 at gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Hi, this is a duplication of
>>>>>>>>
>>>>>>>> https://lists.mcs.anl.gov/pipermail/petsc-users/2022-August/046746.html
>>>>>>>>
>>>>>>>> for in-depth question.
>>>>>>>>
>>>>>>>> I wrote a short code as attached to test interpolation between two
>>>>>>>> DMPlex objects. Goal is to map solution field defined on
>>>>>>>> cell-centroid(DMPlex-1) into vertex(DMPlex-2) field or vice versa.
>>>>>>>>
>>>>>>>> Basically, DMCreateInterpolation() fails in the example code and it
>>>>>>>> was not allowed to get useful error messages. The DMPlex is created by
>>>>>>>> loading a 2D square box in Gmsh. Can I get some comments on that?
>>>>>>>>
>>>>>>>
>>>>>>> Sure, I am looking at it.
>>>>>>>
>>>>>>>   Thanks,
>>>>>>>
>>>>>>>      Matt
>>>>>>>
>>>>>>>
>>>>>>>> Thanks,
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> 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/>
>

-- 
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-users/attachments/20220829/48914cf0/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.F90
Type: application/octet-stream
Size: 10351 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220829/48914cf0/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: x.png
Type: image/png
Size: 194379 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220829/48914cf0/attachment-0001.png>