1ce78bad3SBarry Smith#include <petsc/finclude/petsc.h> 2*c5e229c2SMartin Diehlprogram ex14f90 3ce78bad3SBarry Smith use petsc 4ce78bad3SBarry Smith implicit none 5ce78bad3SBarry Smith 6ce78bad3SBarry Smith type(tDM) :: dm 7ce78bad3SBarry Smith type(tVec) :: u 8ce78bad3SBarry Smith type(tPetscSection) :: section 9ce78bad3SBarry Smith PetscInt :: dim, numFields, numBC 10ce78bad3SBarry Smith PetscMPIInt :: rank 11ce78bad3SBarry Smith PetscInt, dimension(2) :: numComp 12ce78bad3SBarry Smith PetscInt, dimension(12) :: numDof 13ce78bad3SBarry Smith PetscInt, dimension(:), pointer :: remoteOffsets 14ce78bad3SBarry Smith type(tPetscSF) :: pointSF 15ce78bad3SBarry Smith type(tPetscSF) :: sectionSF 16ce78bad3SBarry Smith PetscScalar, dimension(:), pointer :: array 17ce78bad3SBarry Smith PetscReal :: val 18ce78bad3SBarry Smith PetscErrorCode :: ierr 19ce78bad3SBarry Smith PetscInt :: zero = 0, one = 1 20ce78bad3SBarry Smith 21ce78bad3SBarry Smith PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr)) 22ce78bad3SBarry Smith PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr)) 23ce78bad3SBarry Smith PetscCallA(DMSetType(dm, DMPLEX, ierr)) 24ce78bad3SBarry Smith PetscCallA(DMSetFromOptions(dm, ierr)) 25ce78bad3SBarry Smith PetscCallA(DMViewFromOptions(dm, PETSC_NuLL_OBJECT, "-dm_view", ierr)) 26ce78bad3SBarry Smith PetscCallA(DMGetDimension(dm, dim, ierr)) 27ce78bad3SBarry Smith 28ce78bad3SBarry Smith !!! Describe the solution variables that are discretized on the mesh 29ce78bad3SBarry Smith ! Create scalar field u and a vector field v 30ce78bad3SBarry Smith numFields = 2 31ce78bad3SBarry Smith numComp = [one, dim] 32ce78bad3SBarry Smith numDof = 0 33ce78bad3SBarry Smith !Let u be defined on cells 34ccfd86f1SBarry Smith numDof(0*(dim + 1) + dim + 1) = 1 35ce78bad3SBarry Smith !Let v be defined on vertices 36ce78bad3SBarry Smith numDof(1*(dim + 1) + 1) = dim 37ce78bad3SBarry Smith !No boundary conditions */ 38ce78bad3SBarry Smith numBC = 0 39ce78bad3SBarry Smith 40ce78bad3SBarry Smith !!! Create a PetscSection to handle the layout of the discretized variables 41ce78bad3SBarry Smith PetscCallA(DMSetNumFields(dm, numFields, ierr)) 42ce78bad3SBarry Smith 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)) 43ce78bad3SBarry Smith ! Name the Field variables 44ce78bad3SBarry Smith PetscCallA(PetscSectionSetFieldName(section, zero, "u", ierr)) 45ce78bad3SBarry Smith PetscCallA(PetscSectionSetFieldName(section, one, "v", ierr)) 46ce78bad3SBarry Smith ! Tell the DM to use this data layout 47ce78bad3SBarry Smith PetscCallA(DMSetLocalSection(dm, section, ierr)) 48ce78bad3SBarry Smith 49ce78bad3SBarry Smith !!! Construct the communication pattern for halo exchange between local vectors */ 50ce78bad3SBarry Smith ! Get the point SF: an object that says which copies of mesh points (cells, 51ce78bad3SBarry Smith ! vertices, faces, edges) are copies of points on other processes 52ce78bad3SBarry Smith PetscCallA(DMGetPointSF(dm, pointSF, ierr)) 53ce78bad3SBarry Smith ! Relate the locations of ghost degrees of freedom on this process 54ce78bad3SBarry Smith ! to their locations of the non-ghost copies on a different process 55ce78bad3SBarry Smith PetscCallA(PetscSFCreateRemoteOffsets(pointSF, section, section, remoteOffsets, ierr)) 56ce78bad3SBarry Smith ! Use that information to construct a star forest for halo exchange 57ce78bad3SBarry Smith ! for data described by the local section 58ce78bad3SBarry Smith PetscCallA(PetscSFCreateSectionSF(pointSF, section, remoteOffsets, section, sectionSF, ierr)) 59ce78bad3SBarry Smith if (associated(remoteOffsets)) then 60ce78bad3SBarry Smith PetscCallA(PetscSFDestroyRemoteOffsets(remoteOffsets, ierr)) 61ce78bad3SBarry Smith end if 62ce78bad3SBarry Smith 63ce78bad3SBarry Smith !!! Demo of halo exchange 64ce78bad3SBarry Smith ! Create a Vec with this layout 65ce78bad3SBarry Smith PetscCallA(DMCreateLocalVector(dm, u, ierr)) 66ce78bad3SBarry Smith PetscCallA(PetscObjectSetName(u, "Local vector", ierr)) 67ce78bad3SBarry Smith PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr)) 68ce78bad3SBarry Smith ! Set all mesh values to the MPI rank 69ce78bad3SBarry Smith val = rank 70ce78bad3SBarry Smith PetscCallA(VecSet(u, val, ierr)) 71ce78bad3SBarry Smith ! Get the raw array of values 72ce78bad3SBarry Smith PetscCallA(VecGetArrayWrite(u, array, ierr)) 73ce78bad3SBarry Smith !!! HALO EXCHANGE 74ce78bad3SBarry Smith PetscCallA(PetscSFBcastBegin(sectionSF, MPIU_SCALAR, array, array, MPI_REPLACE, ierr)) 75ce78bad3SBarry Smith ! local work can be done between Begin() and End() 76ce78bad3SBarry Smith PetscCallA(PetscSFBcastEnd(sectionSF, MPIU_SCALAR, array, array, MPI_REPLACE, ierr)) 77ce78bad3SBarry Smith ! Restore the raw array of values 78ce78bad3SBarry Smith PetscCallA(VecRestoreArrayWrite(u, array, ierr)) 79ce78bad3SBarry Smith ! View the results: should show which process has the non-ghost copy of each degree of freedom 80ce78bad3SBarry Smith PetscCallA(PetscSectionVecView(section, u, PETSC_VIEWER_STDOUT_WORLD, ierr)) 81ce78bad3SBarry Smith PetscCallA(VecDestroy(u, ierr)) 82ce78bad3SBarry Smith 83ce78bad3SBarry Smith PetscCallA(PetscSFDestroy(sectionSF, ierr)) 84ce78bad3SBarry Smith PetscCallA(PetscSectionDestroy(section, ierr)) 85ce78bad3SBarry Smith PetscCallA(DMDestroy(dm, ierr)) 86ce78bad3SBarry Smith PetscCallA(PetscFinalize(ierr)) 87ce78bad3SBarry Smithend program ex14f90 88ce78bad3SBarry Smith!/*TEST 89ce78bad3SBarry Smith! # Test on a 1D mesh with overlap 90ce78bad3SBarry Smith! test: 91ce78bad3SBarry Smith! nsize: 3 92ce78bad3SBarry Smith! requires: !complex 93ce78bad3SBarry Smith! args: -dm_plex_dim 1 -dm_plex_box_faces 3 -dm_refine_pre 1 -petscpartitioner_type simple -dm_distribute_overlap 1 94ce78bad3SBarry Smith! 95ce78bad3SBarry Smith!TEST*/ 96