[petsc-users] Getting a vector from a DM to output VTK
Nicholas Arnold-Medabalimi
narnoldm at umich.edu
Mon Dec 26 09:48:43 CST 2022
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
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20221226/0d732141/attachment-0001.html>
More information about the petsc-users
mailing list