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