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

[Dharmendar Reddy]

In my code, I apply boundary conditions to the local vector before it
is used for assembling the residual and jacobian.
I use VecSetValuesSection with mode= INSERT_BC_VALUES

I also do the boundary condition's insert on the local vector before
it is used for vecview.

On Wed, Jan 29, 2014 at 12:00 PM, 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