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