1program ex14f90 2 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 7 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 22 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)) 29 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 41 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)) 50 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 64 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)) 84 85 PetscCallA(PetscSFDestroy(sectionSF, ierr)) 86 PetscCallA(PetscSectionDestroy(section, ierr)) 87 PetscCallA(DMDestroy(dm, ierr)) 88 PetscCallA(PetscFinalize(ierr)) 89end 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*/ 101