Actual source code: ex9.c

  1: static char help[] = "Demonstrates use of VecCreateGhost().\n\n";

  3: /*
  4:    Description: Ghost padding is one way to handle local calculations that
  5:       involve values from other processors. VecCreateGhost() provides
  6:       a way to create vectors with extra room at the end of the vector
  7:       array to contain the needed ghost values from other processors,
  8:       vector computations are otherwise unaffected.
  9: */

 11: /*
 12:   Include "petscvec.h" so that we can use vectors.  Note that this file
 13:   automatically includes:
 14:      petscsys.h       - base PETSc routines   petscis.h     - index sets
 15:      petscviewer.h - viewers
 16: */
 17: #include <petscvec.h>

 19: int main(int argc, char **argv)
 20: {
 21:   PetscMPIInt            rank, size;
 22:   PetscInt               nlocal = 6, nghost = 2, ifrom[2], i, rstart, rend;
 23:   PetscBool              flg, flg2, flg3, flg4, flg5;
 24:   PetscScalar            value, *array, *tarray = 0;
 25:   Vec                    lx, gx, gxs;
 26:   IS                     ghost;
 27:   ISLocalToGlobalMapping mapping;

 29:   PetscFunctionBeginUser;
 30:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 31:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 32:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 33:   PetscCheck(size == 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Must run example with two processors");

 35:   /*
 36:      Construct a two dimensional graph connecting nlocal degrees of
 37:      freedom per processor. From this we will generate the global
 38:      indices of needed ghost values

 40:      For simplicity we generate the entire graph on each processor:
 41:      in real application the graph would stored in parallel, but this
 42:      example is only to demonstrate the management of ghost padding
 43:      with VecCreateGhost().

 45:      In this example we consider the vector as representing
 46:      degrees of freedom in a one dimensional grid with periodic
 47:      boundary conditions.

 49:         ----Processor  1---------  ----Processor 2 --------
 50:          0    1   2   3   4    5    6    7   8   9   10   11
 51:                                |----|
 52:          |-------------------------------------------------|

 54:   */

 56:   if (rank == 0) {
 57:     ifrom[0] = 11;
 58:     ifrom[1] = 6;
 59:   } else {
 60:     ifrom[0] = 0;
 61:     ifrom[1] = 5;
 62:   }

 64:   /*
 65:      Create the vector with two slots for ghost points. Note that both
 66:      the local vector (lx) and the global vector (gx) share the same
 67:      array for storing vector values.
 68:   */
 69:   PetscCall(PetscOptionsHasName(NULL, NULL, "-allocate", &flg));
 70:   PetscCall(PetscOptionsHasName(NULL, NULL, "-vecmpisetghost", &flg2));
 71:   PetscCall(PetscOptionsHasName(NULL, NULL, "-minvalues", &flg3));
 72:   if (flg) {
 73:     PetscCall(PetscMalloc1(nlocal + nghost, &tarray));
 74:     PetscCall(VecCreateGhostWithArray(PETSC_COMM_WORLD, nlocal, PETSC_DECIDE, nghost, ifrom, tarray, &gxs));
 75:   } else if (flg2) {
 76:     PetscCall(VecCreate(PETSC_COMM_WORLD, &gxs));
 77:     PetscCall(VecSetType(gxs, VECMPI));
 78:     PetscCall(VecSetSizes(gxs, nlocal, PETSC_DECIDE));
 79:     PetscCall(VecMPISetGhost(gxs, nghost, ifrom));
 80:   } else {
 81:     PetscCall(VecCreateGhost(PETSC_COMM_WORLD, nlocal, PETSC_DECIDE, nghost, ifrom, &gxs));
 82:   }

 84:   /*
 85:       Test VecDuplicate()
 86:   */
 87:   PetscCall(VecDuplicate(gxs, &gx));
 88:   PetscCall(VecDestroy(&gxs));

 90:   /*
 91:      Access the local representation
 92:   */
 93:   PetscCall(VecGhostGetLocalForm(gx, &lx));

 95:   /*
 96:      Set the values from 0 to 12 into the "global" vector
 97:   */
 98:   PetscCall(VecGetOwnershipRange(gx, &rstart, &rend));
 99:   for (i = rstart; i < rend; i++) {
100:     value = (PetscScalar)i;
101:     PetscCall(VecSetValues(gx, 1, &i, &value, INSERT_VALUES));
102:   }
103:   PetscCall(VecAssemblyBegin(gx));
104:   PetscCall(VecAssemblyEnd(gx));

106:   PetscCall(VecGhostUpdateBegin(gx, INSERT_VALUES, SCATTER_FORWARD));
107:   PetscCall(VecGhostUpdateEnd(gx, INSERT_VALUES, SCATTER_FORWARD));

109:   /*
110:      Print out each vector, including the ghost padding region.
111:   */
112:   PetscCall(VecGetArray(lx, &array));
113:   for (i = 0; i < nlocal + nghost; i++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " %g\n", i, (double)PetscRealPart(array[i])));
114:   PetscCall(VecRestoreArray(lx, &array));
115:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
116:   PetscCall(VecGhostRestoreLocalForm(gx, &lx));

118:   /* Another test that sets ghost values and then accumulates onto the owning processors using MIN_VALUES */
119:   if (flg3) {
120:     if (rank == 0) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\nTesting VecGhostUpdate with MIN_VALUES\n"));
121:     PetscCall(VecGhostGetLocalForm(gx, &lx));
122:     PetscCall(VecGetArray(lx, &array));
123:     for (i = 0; i < nghost; i++) array[nlocal + i] = rank ? (PetscScalar)4 : (PetscScalar)8;
124:     PetscCall(VecRestoreArray(lx, &array));
125:     PetscCall(VecGhostRestoreLocalForm(gx, &lx));

127:     PetscCall(VecGhostUpdateBegin(gx, MIN_VALUES, SCATTER_REVERSE));
128:     PetscCall(VecGhostUpdateEnd(gx, MIN_VALUES, SCATTER_REVERSE));

130:     PetscCall(VecGhostGetLocalForm(gx, &lx));
131:     PetscCall(VecGetArray(lx, &array));

133:     for (i = 0; i < nlocal + nghost; i++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " %g\n", i, (double)PetscRealPart(array[i])));
134:     PetscCall(VecRestoreArray(lx, &array));
135:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
136:     PetscCall(VecGhostRestoreLocalForm(gx, &lx));
137:   }

139:   PetscCall(PetscOptionsHasName(NULL, NULL, "-vecghostgetghostis", &flg4));
140:   if (flg4) {
141:     PetscCall(VecGhostGetGhostIS(gx, &ghost));
142:     PetscCall(ISView(ghost, PETSC_VIEWER_STDOUT_WORLD));
143:   }
144:   PetscCall(PetscOptionsHasName(NULL, NULL, "-getgtlmapping", &flg5));
145:   if (flg5) {
146:     PetscCall(VecGetLocalToGlobalMapping(gx, &mapping));
147:     PetscCall(ISLocalToGlobalMappingView(mapping, NULL));
148:   }

150:   PetscCall(VecDestroy(&gx));

152:   if (flg) PetscCall(PetscFree(tarray));
153:   PetscCall(PetscFinalize());
154:   return 0;
155: }

157: /*TEST

159:      test:
160:        nsize: 2

162:      test:
163:        suffix: 2
164:        nsize: 2
165:        args: -allocate
166:        output_file: output/ex9_1.out

168:      test:
169:        suffix: 3
170:        nsize: 2
171:        args: -vecmpisetghost
172:        output_file: output/ex9_1.out

174:      test:
175:        suffix: 4
176:        nsize: 2
177:        args: -minvalues
178:        output_file: output/ex9_2.out
179:        requires: !complex

181:      test:
182:        suffix: 5
183:        nsize: 2
184:        args: -vecghostgetghostis

186:      test:
187:        suffix: 6
188:        nsize: 2
189:        args: -getgtlmapping

191: TEST*/