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

[David Fuentes]

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.

    {



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