[petsc-users] Creating multi-field section using DMPlex

[Shashwat Tiwari]

Hi,
I am working on a Finite Volume scheme for a hyperbolic system of equations
equations with 6 variables on unstructured grids. I am trying to create a
section with 6 fields for this purpose, and have written a small test code
for creating the section by editing the example given at
src/dm/impls/plex/tutorials/ex1.c as follows:

int main(int argc, char **argv)
{
  DM             dm, dmDist = NULL;
  Vec            U;
  PetscSection   section;
  PetscViewer    viewer;
  PetscInt       dim = 2, numFields, numBC, i;
  PetscInt       numComp[6];
  PetscInt       numDof[18];
  PetscBool      interpolate = PETSC_TRUE, useCone = PETSC_TRUE, useClosure
= PETSC_TRUE;
  PetscReal      lower[3], upper[3]; // lower left and upper right
coordinates of domain
  PetscInt       cells[2];
  PetscErrorCode ierr;

  // define domain properties
  cells[0] = 4; cells[1] = 4;
  lower[0] = -1; lower[1] = -1; lower[2] = 0;
  upper[0] =  1; upper[1] =  1; upper[2] = 0;

  ierr = PetscInitialize(&argc, &argv, NULL, help); CHKERRQ(ierr);
  // create the mesh
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Generating mesh ... ");
CHKERRQ(ierr);
  ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE,
                             cells, lower, upper, NULL, interpolate, &dm);
CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Done\n"); CHKERRQ(ierr);

  // set cell adjacency
  ierr = DMSetBasicAdjacency(dm, useCone, useClosure); CHKERRQ(ierr);
  // Distribute mesh over processes
  ierr = DMPlexDistribute(dm, 1, NULL, &dmDist);CHKERRQ(ierr);
  if (dmDist) {ierr = DMDestroy(&dm);CHKERRQ(ierr); dm = dmDist;}

  // Create scalar fields
  numFields  = 6;
  numComp[0] = 1;
  numComp[1] = 1;
  numComp[2] = 1;
  numComp[3] = 1;
  numComp[4] = 1;
  numComp[5] = 1;

  for (i = 0; i < numFields*(dim+1); ++i) numDof[i] = 0;
  numDof[0*(dim+1)+dim] = 1;
  numDof[1*(dim+1)+dim] = 1;
  numDof[2*(dim+1)+dim] = 1;
  numDof[3*(dim+1)+dim] = 1;
  numDof[4*(dim+1)+dim] = 1;
  numDof[5*(dim+1)+dim] = 1;

  numBC = 0;

  // Create a PetscSection with this data layout
  ierr = DMSetNumFields(dm, numFields);CHKERRQ(ierr);
  ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, NULL, NULL,
NULL, NULL, &section);CHKERRQ(ierr);

  // Name the Field variables
  ierr = PetscSectionSetFieldName(section, 0, "h");CHKERRQ(ierr);
  ierr = PetscSectionSetFieldName(section, 1, "m1");CHKERRQ(ierr);
  ierr = PetscSectionSetFieldName(section, 2, "m2");CHKERRQ(ierr);
  ierr = PetscSectionSetFieldName(section, 3, "E11");CHKERRQ(ierr);
  ierr = PetscSectionSetFieldName(section, 4, "E12");CHKERRQ(ierr);
  ierr = PetscSectionSetFieldName(section, 5, "E22");CHKERRQ(ierr);

  // set adjacency for each field
  ierr = DMSetAdjacency(dm, 0, useCone, useClosure); CHKERRQ(ierr);
  ierr = DMSetAdjacency(dm, 1, useCone, useClosure); CHKERRQ(ierr);
  ierr = DMSetAdjacency(dm, 2, useCone, useClosure); CHKERRQ(ierr);
  ierr = DMSetAdjacency(dm, 3, useCone, useClosure); CHKERRQ(ierr);
  ierr = DMSetAdjacency(dm, 4, useCone, useClosure); CHKERRQ(ierr);
  ierr = DMSetAdjacency(dm, 5, useCone, useClosure); CHKERRQ(ierr);

  ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  // Tell the DM to use this data layout
  ierr = DMSetLocalSection(dm, section);CHKERRQ(ierr);
  ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  // Create a Vec with this layout and view it
  ierr = DMGetGlobalVector(dm, &U);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) U, "U"); CHKERRQ(ierr);
  ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);
  ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
  ierr = PetscViewerPushFormat(viewer,
PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
  ierr = PetscViewerFileSetName(viewer, "sol.vtk");CHKERRQ(ierr);
  //ierr = PetscViewerPushFormat(viewer,
PETSC_VIEWER_VTK_VTU);CHKERRQ(ierr);
  //ierr = PetscViewerFileSetName(viewer, "sol.vtu");CHKERRQ(ierr);
  ierr = VecView(U, viewer);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(dm, &U);CHKERRQ(ierr);

  // Cleanup
  ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

When exporting the vector data to "vtk" format, the code works fine and I
am able to see the 6 field names in the visualizer. But when I change it to
"vtu", I get the following error:

[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Total number of field components 1 != block size 6
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.13.0, unknown
[0]PETSC ERROR: ./ex1 on a arch-linux2-c-debug named shashwat by shashwat
Fri May 22 22:11:13 2020
[0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++
--with-fc=gfortran --download-mpich --download-fblaslapack
--with-debugging=1 --download-triangle --download-metis --download-parmetis
--download-cmake
[0]PETSC ERROR: #1 DMPlexVTKWriteAll_VTU() line 334 in
/home/shashwat/local/petsc/src/dm/impls/plex/plexvtu.c
[0]PETSC ERROR: #2 DMPlexVTKWriteAll() line 688 in
/home/shashwat/local/petsc/src/dm/impls/plex/plexvtk.c
[0]PETSC ERROR: #3 PetscViewerFlush_VTK() line 100 in
/home/shashwat/local/petsc/src/sys/classes/viewer/impls/vtk/vtkv.c
[0]PETSC ERROR: #4 PetscViewerFlush() line 26 in
/home/shashwat/local/petsc/src/sys/classes/viewer/interface/flush.c
[0]PETSC ERROR: #5 PetscViewerDestroy() line 113 in
/home/shashwat/local/petsc/src/sys/classes/viewer/interface/view.c
[0]PETSC ERROR: #6 main() line 485 in ex1.c
[0]PETSC ERROR: No PETSc Option Table entries
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to petsc-maint at mcs.anl.gov----------
application called MPI_Abort(MPI_COMM_SELF, 485077) - process 0

I would like to use the "vtu" format for visualiztion. Kindly have a look
and let me know I might be doing wrong here, and how can I make it work.

Regards,
Shashwat
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200522/7d3b741b/attachment.html>