[petsc-users] Getting a vector from a DM to output VTK

[Matthew Knepley]

On Mon, Dec 26, 2022 at 3:21 AM Nicholas Arnold-Medabalimi <
narnoldm at umich.edu> wrote:

> Hi Matt
>
> I was able to get this all squared away. It turns out I was initializing
> the viewer incorrectly—my mistake. However, there is a follow-up question.
> A while back, we discussed distributing a vector field from an initial DM
> to a new distributed DM. The way you said to do this was
>
>     // Distribute the submesh with overlap of 1
>     DMPlexDistribute(sub_da, overlap, &distributionSF, &sub_da_dist);
>     //Create a vector and section for the distribution
>     VecCreate(PETSC_COMM_WORLD, &state_dist);
>     VecSetDM(state_dist, sub_da_dist);
>     PetscSectionCreate(PETSC_COMM_WORLD, &distSection);
>     DMSetLocalSection(sub_da_dist, distSection);
>     DMPlexDistributeField(sub_da_dist, distributionSF, filteredSection,
> state_filtered, distSection, state_dist);
>   I've forgone Fortran to debug this all in C and then integrate the
> function calls into the Fortran code.
>
> There are two questions here.
>
> 1) How do I associate a vector associated with a DM using VecSetDM to
> output properly as a VTK? When I call VecView at present, if I call VecView
> on state_dist, it will not output anything.
>

This is a problem. The different pieces of interface were added at
different times. We should really move that manipulation of the function
table into VecSetDM(). Here is the code:


https://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/plexcreate.c#L4135

You can make the overload call yourself for now, until we decide on the
best fix.


> 2) The visualization is nice, but when I look at the Vec of the
> distributed field using stdout, something isn't distributing correctly, as
> the vector still has some uninitialized values. This is apparent if I
> output the original vector and the distributed vector. Examining the inside
> of DMPlexDistributeField I suspect I'm making a mistake with the sections
> I'm passing. filtered section in this case is the global section but if I
> try the local section I get an error so I'm not sure.
>

These should definitely be local sections. Global sections are always built
after the fact, and building the global section needs the SF that indicates
what points are shared, not the distribution SF that moves points. I need
to go back and put in checks that all the arguments are the right type.
Thanks for bringing that up. Lets track down the error for local sections.

  Matt


> *Original Vector(state_filtered)*
> Vec Object: Vec_0x84000004_1 2 MPI processes
>   type: mpi
> Process [0]
> 101325.
> 300.
> 101326.
> 301.
> 101341.
> 316.
> Process [1]
> 101325.
> 300.
> 101326.
> 301.
> 101345.
> 320.
> 101497.
> 472.
> 101516.
> 491.
> *Re-Distributed Vector (state_dist) *
> Vec Object: 2 MPI processes
>   type: mpi
> Process [0]
> 101325.
> 300.
> 101326.
> 301.
> 101341.
> 316.
> 7.90505e-323
> 1.97626e-323
> 4.30765e-312
> 6.91179e-310
> Process [1]
> 101497.
> 472.
> 101516.
> 491.
> 1.99665e-314
> 8.14714e-321
>
>
> Any insight on distributing this field and correcting the error would be
> appreciated.
>
> Sincerely and Happy Holiday
> Nicholas
>
>
>
> On Fri, Dec 23, 2022 at 10:57 AM Matthew Knepley <knepley at gmail.com>
> wrote:
>
>> On Thu, Dec 22, 2022 at 12:41 AM Nicholas Arnold-Medabalimi <
>> narnoldm at umich.edu> wrote:
>>
>>> Hi Petsc Users
>>>
>>> I've been having trouble consistently getting a vector generated from a
>>> DM to output to VTK correctly. I've used ex1.c (which works properly)to try
>>> and figure it out, but I'm still having some issues. I must be missing
>>> something small that isn't correctly associating the section with the DM.
>>>
>>>     DMPlexGetChart(dm, &p0, &p1);
>>>     PetscSection section_full;
>>>     PetscSectionCreate(PETSC_COMM_WORLD, &section_full);
>>>     PetscSectionSetNumFields(section_full, 1);
>>>     PetscSectionSetChart(section_full, p0, p1);
>>>     PetscSectionSetFieldName(section_full, 0, "state");
>>>
>>>     for (int i = c0; i < c1; i++)
>>>     {
>>>         PetscSectionSetDof(section_full, i, 1);
>>>         PetscSectionSetFieldDof(section_full, i, 0, 1);
>>>     }
>>>     PetscSectionSetUp(section_full);
>>>     DMSetNumFields(dm, 1);
>>>     DMSetLocalSection(dm, section_full);
>>>     DMCreateGlobalVector(dm, &state_full);
>>>
>>>     int o0, o1;
>>>     VecGetOwnershipRange(state_full, &o0, &o1);
>>>     PetscScalar *state_full_array;
>>>     VecGetArray(state_full, &state_full_array);
>>>
>>>     for (int i = 0; i < (c1 - c0); i++)
>>>     {
>>>         int offset;
>>>         PetscSectionGetOffset(section_full, i, &offset);
>>>         state_full_array[offset] = 101325 + i;
>>>     }
>>>
>>>     VecRestoreArray(state_full, &state_full_array);
>>>
>>>
>>>     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
>>>     PetscViewerSetType(viewer, PETSCVIEWERVTK);
>>>     PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
>>>     PetscViewerFileSetName(viewer, "mesh.vtu");
>>>     VecView(state_full, viewer);
>>>
>>> If I run this mesh.vtu isn't generated at all. If I instead do a DMView
>>> passing the DM it will just output the mesh correctly.
>>>
>>> Any assistance would be greatly appreciated.
>>>
>>
>> DMCreateGlobalVector() dispatches to DMCreateGlobalVector_Plex(), which
>> resets the view method to VecView_Plex(), which should dispatch to
>> VecView_Plex_Local_VTK(). You can verify this in the debugger, or send us
>> code we can run to verify it.
>>
>>   Thanks,
>>
>>     Matt
>>
>>
>>> Sincerely
>>> Nicholas
>>>
>>>
>>> --
>>> Nicholas Arnold-Medabalimi
>>>
>>> Ph.D. Candidate
>>> Computational Aeroscience Lab
>>> University of Michigan
>>>
>>
>>
>> --
>> 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/>
>>
>
>
> --
> Nicholas Arnold-Medabalimi
>
> Ph.D. Candidate
> Computational Aeroscience Lab
> University of Michigan
>


-- 
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/20221226/3664a572/attachment-0001.html>