Actual source code: ex14f90.F90

      1: program ex14f90

      3: #include <petsc/finclude/petsc.h>
      4:   use petsc
      5:   use mpi     ! needed when PETSC_HAVE_MPI_F90MODULE is not true to define MPI_REPLACE
      6:   implicit none

      8:   type(tDM)                        :: dm
      9:   type(tVec)                       :: u
     10:   type(tPetscSection)              :: section
     11:   PetscInt                         :: dim, numFields, numBC
     12:   PetscMPIInt                      :: rank
     13:   PetscInt, dimension(2)            :: numComp
     14:   PetscInt, dimension(12)           :: numDof
     15:   PetscInt, dimension(:), pointer    :: remoteOffsets
     16:   type(tPetscSF)                   :: pointSF
     17:   type(tPetscSF)                   :: sectionSF
     18:   PetscScalar, dimension(:), pointer :: array
     19:   PetscReal                        :: val
     20:   PetscErrorCode                   :: ierr
     21:   PetscInt                         :: zero = 0, one = 1

     23:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
     24:   PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr))
     25:   PetscCallA(DMSetType(dm, DMPLEX, ierr))
     26:   PetscCallA(DMSetFromOptions(dm, ierr))
     27:   PetscCallA(DMViewFromOptions(dm, PETSC_NuLL_OBJECT, "-dm_view", ierr))
     28:   PetscCallA(DMGetDimension(dm, dim, ierr))

     30:   !!! Describe the solution variables that are discretized on the mesh
     31:   ! Create scalar field u and a vector field v
     32:   numFields = 2
     33:   numComp = [one, dim]
     34:   numDof = 0
     35:   !Let u be defined on cells
     36:   numDof(0*(dim + 1) + dim + 1) = 1
     37:   !Let v be defined on vertices
     38:   numDof(1*(dim + 1) + 1) = dim
     39:   !No boundary conditions */
     40:   numBC = 0

     42:   !!! Create a PetscSection to handle the layout of the discretized variables
     43:   PetscCallA(DMSetNumFields(dm, numFields, ierr))
     44:   PetscCallA(DMPlexCreateSection(dm, PETSC_NULL_DMLABEL_ARRAY, numComp, numDof, numBC, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_IS_ARRAY, PETSC_NULL_IS_ARRAY, PETSC_NULL_IS, section, ierr))
     45:   ! Name the Field variables
     46:   PetscCallA(PetscSectionSetFieldName(section, zero, "u", ierr))
     47:   PetscCallA(PetscSectionSetFieldName(section, one, "v", ierr))
     48:   ! Tell the DM to use this data layout
     49:   PetscCallA(DMSetLocalSection(dm, section, ierr))

     51:   !!! Construct the communication pattern for halo exchange between local vectors */
     52:   ! Get the point SF: an object that says which copies of mesh points (cells,
     53:   ! vertices, faces, edges) are copies of points on other processes
     54:   PetscCallA(DMGetPointSF(dm, pointSF, ierr))
     55:   ! Relate the locations of ghost degrees of freedom on this process
     56:   ! to their locations of the non-ghost copies on a different process
     57:   PetscCallA(PetscSFCreateRemoteOffsets(pointSF, section, section, remoteOffsets, ierr))
     58:   ! Use that information to construct a star forest for halo exchange
     59:   ! for data described by the local section
     60:   PetscCallA(PetscSFCreateSectionSF(pointSF, section, remoteOffsets, section, sectionSF, ierr))
     61:   if (associated(remoteOffsets)) then
     62:     PetscCallA(PetscSFDestroyRemoteOffsets(remoteOffsets, ierr))
     63:   end if

     65:   !!! Demo of halo exchange
     66:   ! Create a Vec with this layout
     67:   PetscCallA(DMCreateLocalVector(dm, u, ierr))
     68:   PetscCallA(PetscObjectSetName(u, "Local vector", ierr))
     69:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     70:   ! Set all mesh values to the MPI rank
     71:   val = rank
     72:   PetscCallA(VecSet(u, val, ierr))
     73:   ! Get the raw array of values
     74:   PetscCallA(VecGetArrayWrite(u, array, ierr))
     75:   !!! HALO EXCHANGE
     76:   PetscCallA(PetscSFBcastBegin(sectionSF, MPIU_SCALAR, array, array, MPI_REPLACE, ierr))
     77:   ! local work can be done between Begin() and End()
     78:   PetscCallA(PetscSFBcastEnd(sectionSF, MPIU_SCALAR, array, array, MPI_REPLACE, ierr))
     79:   ! Restore the raw array of values
     80:   PetscCallA(VecRestoreArrayWrite(u, array, ierr))
     81:   ! View the results: should show which process has the non-ghost copy of each degree of freedom
     82:   PetscCallA(PetscSectionVecView(section, u, PETSC_VIEWER_STDOUT_WORLD, ierr))
     83:   PetscCallA(VecDestroy(u, ierr))

     85:   PetscCallA(PetscSFDestroy(sectionSF, ierr))
     86:   PetscCallA(PetscSectionDestroy(section, ierr))
     87:   PetscCallA(DMDestroy(dm, ierr))
     88:   PetscCallA(PetscFinalize(ierr))
     89: end program ex14f90
     90: !/*TEST
     91: !  build:
     92: !    requires: defined(PETSC_USING_F90FREEFORM)
     93: !
     94: !  # Test on a 1D mesh with overlap
     95: !  test:
     96: !    nsize: 3
     97: !    requires: !complex
     98: !    args: -dm_plex_dim 1 -dm_plex_box_faces 3 -dm_refine_pre 1 -petscpartitioner_type simple -dm_distribute_overlap 1
     99: !
    100: !TEST*/