[petsc-users] Creating multi-field section using DMPlex
Matthew Knepley
knepley at gmail.com
Fri May 22 20:07:27 CDT 2020
On Fri, May 22, 2020 at 12:57 PM Shashwat Tiwari <shaswat121994 at gmail.com>
wrote:
> 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:
>
If you use the PetscFV stuff, it assumes everything is in a single field.
This is because FV methods often want Riemann solvers with all the fields,
but Plex
splits the callbacks by field in order to get sparsity in the finite
element matrix. In order to make things more understandable, we added
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PetscSection/PetscSectionSetComponentName.html#PetscSectionSetComponentName
so you can name the components just like you would name fields. I think
this may be what you want. Let me know if this does not help.
Thanks,
Matt
> 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, §ion);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(§ion);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
>
--
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
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200522/8e44859c/attachment.html>
More information about the petsc-users
mailing list