[petsc-users] Getting a vector from a DM to output VTK
Matthew Knepley
knepley at gmail.com
Mon Dec 26 09:51:51 CST 2022
On Mon, Dec 26, 2022 at 10:49 AM Nicholas Arnold-Medabalimi <
narnoldm at umich.edu> wrote:
> Oh it's not a worry. I'm debugging this first in C++, and once it's
> working I don't actually need to view what's happening in Fortran when I
> move over. In my C debugging code. After I create the distribution vector
> and distribute the field based on your input, I'm adding
>
> VecSetOperation(state_dist, VECOP_VIEW, (void(*)(void))VecView_Plex);
>
> But I can't compile it because VecView_Plex is undefined. Thanks
>
You need
#include <petsc/private/dmpleximpl.h>
Thanks
Matt
> Sincerely
> Nicholas
>
>
>
> On Mon, Dec 26, 2022 at 10:45 AM Matthew Knepley <knepley at gmail.com>
> wrote:
>
>> On Mon, Dec 26, 2022 at 10:40 AM Nicholas Arnold-Medabalimi <
>> narnoldm at umich.edu> wrote:
>>
>>> Hi Matt
>>>
>>> 1) I'm not sure I follow how to call this. If I insert the
>>> VecSetOperation call I'm not exactly sure what the VecView_Plex is or where
>>> it is defined?
>>>
>>
>> Shoot, this cannot be done in Fortran. I will rewrite this step to get it
>> working for you. It should have been done anyway. I cannot do it
>> today since I have to grade finals, but I should have it this week.
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> 2) Otherwise I've solved this problem with the insight you provided into
>>> the local section. Things look good on the ASCII output but if we can
>>> resolve 1 then I think the loop is fully closed and I can just worry about
>>> the fortran translation.
>>>
>>> Thanks again for all your help.
>>>
>>> Sincerely
>>> Nicholas
>>>
>>>
>>>
>>> On Mon, Dec 26, 2022 at 9:37 AM Matthew Knepley <knepley at gmail.com>
>>> wrote:
>>>
>>>> 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, §ion_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/>
>>>>
>>>
>>>
>>> --
>>> 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/978c89f2/attachment-0001.html>
More information about the petsc-users
mailing list