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

[Matthew Knepley]

On Mon, Mar 31, 2014 at 9:27 AM, David Fuentes <fuentesdt at gmail.com> wrote:

> Hi,
>
> I'm trying to get DMPlex to solve and display with Dirichlet boundary
> conditions.
>
> My global and local solution vectors are of different sizes. The length
> difference is the number of nodes read from the nodeset read in the exodus
> file.
> I am able to project the initial solution/boundary data that I want with
> with VecSetValuesSection. The initial condition looks ok with the global
> solution
> in the viewer. However, I seem to be losing the dirichlet boundary data on
> the solve ? I'm trying to use the suggested global to local scattering on
> the output.
>

I have just worked all this out in 'next', and most of it is in 'master'
now. I will talk through the process

1) Suppose you view a global vector, either VTK or HDF5

2) PETSc maps this to a local vector, and calls the local version, which is
the same as viewing a local vector

3) PETSc will insert the correct boundary values (via
DMPlexInsertBoundaryValuesFEM) if you specify them
    using the Section interface + DMAddBoundary(). Otherwise, you can view
a local vector with the correct boundary values.

4) A scatter is created to a global vector including the boundary values

5) This expanded global vector is passes to the normal HDF5/VTK routines

It is still a little rough, but it should not be much effort to get it
working the way you want.

   Matt


>     {
>
>
>
>  Vec            localX;
>
>
>
> PetscSection   section;
>
>
>
>  PetscInt       vStart, vEnd, v ;
>
>       DMGetLocalVector(dm, &localX);
>       DMGetDefaultSection(dm, &section);
>
>       DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
>       for (v = vStart; v < vEnd; ++v) {
>         PetscScalar bcData = 101.0;
>         VecSetValuesSection(localX, section, v, &bcData ,
> INSERT_BC_VALUES);
>       }
>       DMRestoreLocalVector(dm, &localX);
>       DMLocalToGlobalBegin(dm, localX, INSERT_VALUES, u);
>       DMLocalToGlobalEnd(dm, localX, INSERT_VALUES, u);
>
>       PetscViewer viewer;
>       char               vtkfilename[PETSC_MAX_PATH_LEN] =
> "initialcondition.vtk";
>       ierr =
> PetscViewerVTKOpen(PETSC_COMM_WORLD,vtkfilename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
>       ierr = VecView(u,viewer);CHKERRQ(ierr);
>       ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr)
>     }
>
>   .
>   .
>   ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);
>   .
>   .
>
>   if (write_output) {
>     PetscViewer viewer;
>     const char *name;
>     Vec            lv;
>     DMGetLocalVector(dm, &lv);
>     ierr = PetscObjectGetName((PetscObject)  u,&name);CHKERRQ(ierr);
>     ierr = PetscObjectSetName((PetscObject) lv, name);CHKERRQ(ierr);
>     ierr = DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);CHKERRQ(ierr);
>     ierr = DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);CHKERRQ(ierr);
>     char               vtkfilename[PETSC_MAX_PATH_LEN] = "solution.vtk";
>     ierr =
> PetscViewerVTKOpen(PETSC_COMM_WORLD,vtkfilename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
>     ierr = VecView(lv,viewer);CHKERRQ(ierr);
>     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
>     ierr = DMRestoreLocalVector(dm, &lv);CHKERRQ(ierr);
>   }
>
>
>
> On Wed, Jan 29, 2014 at 2:04 PM, Geoffrey Irving <irving at naml.us> wrote:
>
>> On Wed, Jan 29, 2014 at 12:00 PM, Jed Brown <jed at jedbrown.org> wrote:
>> > Local vectors are supposed to be just that: local. VecViewing a local
>> vector and expecting it to be parallel is perverse.  So we need a real
>> interface.
>> >
>> > Placing 1.0 on the diagonal (and don't assemble into those rows and
>> columns) is the common way to deal with Dirichlet boundary nodes. See ex48
>> for one example. I have written about this in a few places; I can find the
>> more complete description when I have a keyboard.
>>
>> I'll look forward to the improved interface.  For better or worse, I'd
>> like to be able to view simulations in the near future, so it looks
>> like I have to go with Matt's perverse but working version now.
>>
>> Geoffrey
>>
>> > On January 29, 2014 12:53:26 PM MST, Geoffrey Irving <irving at naml.us>
>> wrote:
>> >>On Wed, Jan 29, 2014 at 11:36 AM, Jed Brown <jed at jedbrown.org> wrote:
>> >>> 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).
>> >>
>> >>Matt's sample code doesn't set it either, but for Matt's sample code I
>> >>know where to insert the one line call to DMPlexProjectFunctionLocal.
>> >>For your version I never have explicit access to the local vector, so
>> >>I can't insert the fix.
>> >>
>> >>> 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.
>> >>
>> >>To leave the boundary conditions in, we would need efficient support
>> >>for a very large, very sparse MatNullSpace.  This is doable via
>> >>shells, but is it easy to do in a way that doesn't interfere with the
>> >>user's other null spaces?
>> >>
>> >>Geoffrey
>> >
>>
>
>


-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20140331/4bfbc392/attachment.html>