[petsc-dev] getting data out of a PetscFE simulation

[Geoffrey Irving]

On Tue, Jan 28, 2014 at 8:57 PM, Geoffrey Irving <irving at naml.us> wrote:
> On Tue, Jan 28, 2014 at 8:24 PM, Jed Brown <jed at jedbrown.org> wrote:
>> Matthew Knepley <knepley at gmail.com> writes:
>>
>>> On Tue, Jan 28, 2014 at 9:41 PM, Jed Brown <jed at jedbrown.org> wrote:
>>>
>>>> Matthew Knepley <knepley at gmail.com> writes:
>>>> > In ex62:
>>>> >
>>>> >   if (user.runType == RUN_FULL) {
>>>> >     PetscViewer viewer;
>>>> >     Vec         uLocal;
>>>> >     const char *name;
>>>> >
>>>> >     ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);
>>>> >     ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
>>>> >     ierr = PetscViewerSetFormat(viewer,
>>>> > PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
>>>> >     ierr = PetscViewerFileSetName(viewer, "ex62_sol.vtk");CHKERRQ(ierr);
>>>>
>>>> Please don't do this.  Just set the type to PETSCVIEWERVTK, set the file
>>>> name to some *.vtu; the VTK viewer will infer output format from
>>>> extension and write the VTU XML format with binary-appended data in a
>>>> fairly fast and memory-scalable way.  Then
>>>>
>>>>   VecView(u,viewer);
>>>
>>> I think Jed is wrong here. You need to be using local vectors since
>>> they have BC and constraints in them. Please explain how your VTU
>>> format will magically put them in.
>>
>> We made this work when I threw a temper tantrum while we were writing TS
>> ex11.c (late 2012).
>>
>> PetscErrorCode VecView_Plex(Vec v, PetscViewer viewer)
>> {
>>   DM             dm;
>>   PetscBool      isvtk;
>>   PetscErrorCode ierr;
>>
>>   PetscFunctionBegin;
>>   ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
>>   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
>>   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
>>   if (isvtk) {
>>     Vec         locv;
>>     const char *name;
>>
>>     ierr = DMGetLocalVector(dm, &locv);CHKERRQ(ierr);
>>     ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr);
>>     ierr = PetscObjectSetName((PetscObject) locv, name);CHKERRQ(ierr);
>>     ierr = DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr);
>>     ierr = DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr);
>>     ierr = VecView_Plex_Local(locv, viewer);CHKERRQ(ierr);
>>     ierr = DMRestoreLocalVector(dm, &locv);CHKERRQ(ierr);
>
> Which of those lines adds the correct boundary condition values?

It looks like one of the following is true:

1. VecView on a DMPlex global vector does not insert correct boundary
conditions.
2. Scalar, residual, and jacobian integration are all performing an
extra boundary condition enforcement step.

I'm guessing it's (1).  If so, does that mean I have to do it Matt's
more verbose way?

Geoffrey