[petsc-users] Getting a vector from a DM to output VTK
Matthew Knepley
knepley at gmail.com
Mon Dec 26 09:44:55 CST 2022
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20221226/6a96b922/attachment.html>
More information about the petsc-users
mailing list