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

[Jed Brown]

Matt's sample code doesn't set it either. We need to fix (and the *only* acceptable fix if that VecView always does the right thing, because we have to be able to call it in analysis settings that know nothing about your discretization).

The problem is that some vectors reside in a homogeneous space (e.g. increments and eigenvectors) while others reside in the inhomogeneous space (solutions). We can add a flag or BC attribute on the vector to this effect, but this (and slip conditions) was the issue that led me to conclude that removing boundary nodes was mostly a false economy.

On January 29, 2014 11:00:51 AM MST, Geoffrey Irving <irving at naml.us> wrote:
>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