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

[Nicholas Arnold-Medabalimi]

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.

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.


*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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20221226/937669ba/attachment.html>